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SECTION  1 


INTRODUCTION 

The  increasing  use  of  composite  materials  in  structural  applications,  such  as 

automobiles,  aircraft  and  space  structures,  is  characterized  by  their  high  strength 
(stiH'ness)-to-w eight  ratio,  low  maintenance  costs  and  the  flexibility  in  tailoring  the 
stiffness  and  strength  to  design  requirements.  As  fiber  reinforced  laminates  have 

played  a  more  important  role  in  high  performance  structures  for  the  last  2  decades,  the 
need  to  have  accurate  stress  and  failure  analysis  become  apparent  for  design  or  repair 
purpose. 

Recent  development  in  the  analysis  of  composite  laminate  coupons  under  uniform 
extension  indicated  that  the  high  interlaminar  stresses  near  the  free-edge  are  mainly 
responsible  for  delamination  failure  [l ].  Before  delamination  can  be  predicted  on  the 
basis  of  a  stress-based  failure  criterion,  it  is  essential  that  a  highly  reliable  estimate  of 
interlaminar  stresses  be  available  for  the  given  situation.  However,  it  has  been 
difficult  to  obtain  solutions  for  the  stress  field  Irecause  of  ’he  anisotropy  as  well  as 
heterogeneity  of  the  material,  and  the  difficulty  of  satisfying  traction-free  boundary 

condition  in  a  solution  procedure  based  on  the  displacement  formulation. 

Considerable  research  efforts  have  been  devoted  to  the  study  of  such  free-edge 

delamination  problem.  These  can  be  classified  as  analytical  and  numerical  approaches. 
The  analytical  solutions  are  based  upon  simple  elastic  approximation  [2,3],  modified 
higher  order  theory  [4],  Galerkin  method  [5],  Perturbation  technique  [6],  Boundary  layer 
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theory  [7],  Reissner's  variational  principle  [8,9],  Global-local  model  [  1  ()]  etc.,  while  the 
numerical  solutions  are  based  on  finite  difference  [  1 1,12]  and  finite  element  methods 
including  displacement  [13-16],  stress  [17]  and  hybrid  [18,19]  formulations.  It  was 
found  that  some  of  the  solution  techniques  were  only  applicable  under  certain 
conditions.  For  this  reason,  a  complete  stress  distribution  was  usually  hard  to  obtain. 
Although  results  calculated  from  various  approaches  have  demonstrated  similarities  in 
some  cases,  discrepancies  do  exist  in  the  magnitude  as  well  as  sign  of  the  computed 
interlaminar  stresses  near  the  free-edge  of  laminate  coupons.  One  example  is  shown  in 
Figure  (1)  in  which  significant  difference  was  observed  for  cr,  stress  distribution  along 
the  interface  of  [45/-45],  laminate  based  on  various  solution  techniques  [19], 
Apparently,  one  possible  source  of  these  discrepancies  is  that,  in  these  methods,  the 
continuity  conditions  for  displacements  and  tractions  across  laminate  interfaces  along 
with  traction-free  boundary  condition  along  free-edges  characteristic  of  the  real  life 
situation,  can  only  be  approximated  to  a  limited  extent.  However,  the  credibility  of 
various  methods  in  predicting  the  cr^  distribution  shown  in  Figure  (l)  will  be  judged 
later. 

Due  to  the  presence  of  singular  interlaminar  stresses  near  the  laminate 
free-boundary,  edge  delamination  associated  with  various  types  of  damage  modes,  such 
as  fiber  breakage,  matrix  cracking,  fiber-matrix  debonding,  etc.,  are  observed  to  occur 
under  incremental  loading.  Delamination  can  be  simply  interpreted  as  separation  of 
laminae  from  each  other  in  the  laminate,  and  can  occur  under  static,  impact  or  fatigue 
loading  conditions.  For  the  case  of  a  symmetric  laminate  under  inplane  loading,  the 
strain  components  are  essentially  uniform  throughout  the  laminate.  Due  to  the 
free-edge  effect  the  out-of-plane  interlaminar  stresses,  however,  may  be  sufficiently 
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Figure  1: 


Comparison  of  Z-stress  at  the  Interface  of  45/-45  Layers 


large  to  damage  the  matrix  material,  which  bonds  adjacent  plies  together,  and  cause 
delamination. 

Two  generic  approaches  are  available  for  investigating  damage  modes  in  composite 
materials.  The  first  approach  involves  a  detailed  stress  analysis  used  in  conjunction 
with  a  failure  criterion  to  predict,  and  measure  experimentally,  the  onset  of  fiber 
fracture,  matrix  cracking  and  delamination.  This  can  be  referred  to  as  the  strength 
characterization  approach.  In  the  second  approach,  classical  linear  elastic  fracture 
mechanics  can  be  applied  to  characterize  matrix  cracking  and  the  delamination  process. 
Delamination  has  usually  been  isolated  from  the  other  damage  modes  and  treated  as  a 
stable  crack  growth  [20-25],  and  the  basic  character  of  the  strain  energy  release  rate 
has  been  widely  used  to  predict  the  kinematic  behavior  of  delamination.  However, 
experiments  [26]  have  indicated  that  delamination  usually  does  not  produce  a  clean 
surface  between  the  adjacent  plies;  instead  it  is  associated  with  other  types  of  damage 
such  as  matrix  cracking  and  fiber  breakage.  Thus,  use  of  linear  elastic  fracture 
mechanics  approach  to  study  delamination  growth  seems  to  be  inappropriate. 
Meanwhile,  due  to  the  irregular  occurrence  of  various  damage  modes  in  the  form  of 
different  cracking  patterns,  use  of  anisotropic  strength  and  failure  criteria  is  apparently 
superior  to  the  fracture  mechanic  approach  for  the  determination  of  damage 
characteristics  such  as  type  of  failure  mode,  damage  zone  and  crack  growth  behavior 
including  delamination. 

The  primary  objective  of  the  present  research  was  to  develop  a  finite  element 
model  with  a  sound  theoretical  background,  which  could  accurately  and  efficiently 
predict  the  complete  stress  field  of  the  free-edge  stress  problem  in  composite  laminates 
without  resorting  to  any  special  singularity  elements.  The  next  was  to  incorporate 
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various  commonly  known  macroscopic  failure  criteria  into  the  finite  element 
computational  procedure  to  evaluate  the  performance  of  various  criteria  on  the 
determination  of  onset  of  matrix  cracking  and  delamination  in  the  composite  laminate 
specimens  under  uniform  extension.  In  Section  II,  a  review  of  analytical  and 
numerical  methods  related  to  the  free-edge  stress  problem  is  presented.  Section  111 
contains  the  theoretical  foundation  of  the  finite  element  formulation  including  basic 
variational  principles.  Section  IV  describes  a  continuous  strain  finite  element  model 
based  on  a  compatible  cubic  interpolation  function.  A  continuous  traction  finite 
element  procedure  for  analysis  of  free-edge  delamination  specimens  is  developed  in 
Section  V.  In  Section  VI,  analysis  of  free-edge  effect  as  well  as  onset  of  delamination 
in  various  types  of  laminated  specimens  are  presented.  Section  VII  contains  discussion 
of  the  proposed  finite  element  models  for  analysis  of  free-edge  delamination  specimens. 
Derivation  of  Felippa's  compatible  cubic  interpolation  function  is  summarized  in  the 
Appendix. 
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SECTION  II 


REVIEW  OF  EARLIER  WORK 

2.1  Introduction 

The  problem  of  calculating  interlaminar  '-tresses  near  the  tree-edges  of  a  lave  red 

composite  under  uniform  inplane  extension  has  ireen  investigated  by  many  researchers. 
Most  approximate  solutions  [2-7,10-19]  are  based  upon  elasticity  theory  and  treat  the 
problem  as  a  generalized  plane  strain  case.  This  is  because  first  of  all,  the  classical 

and  even  many  of  the  refined  laminate  theories,  are  single-layer  theories  which  do  not 
account  for  local  effects  such  as  geometric  and  material  discontinuities,  and  the  presence 
of  a  free-edge;  secondly,  use  of  discrete  layered  theory  is  very  uneconomical  and 
impractical  from  the  computational  stand  point.  An  effective  modulus  formulation  [27] 
in  which  each  layer  is  characterized  as  a  homogeneous,  anisotropic  material  has  been 
widely  used  [1-19].  A  complex  state  of  stress  with  high  gradients  has  been  noticed 

[19]  in  the  neighborhood  of  the  free-edge  due  to  the  presence  of  interlaminar  stresses  to 
keep  the  laminae  in  a  state  of  equilibrium.  In  order  to  have  a  precise  prediction  of 

delamination  behavior,  an  accurate  estimate  for  the  near-field  stress  distribution  is 
essential.  However,  due  to  the  singular  nature  of  the  boundary-layer  stress  field  [19], 
an  exact  solution  is  currently  unavailable,  and  discrepancies  exist  in  the  magnitude  and 
even  the  sign  of  the  computed  interlaminar  stresses  near  the  free-edge  (Figure  1)  based 
on  various  approximate  theories. 
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2.2  Analytical  Approach 


Except  for  Pagano's  [8]  approximate  theory  based  on  Reissner's  variational  principle 
and  Pagano  and  Soni's  [l 0]  Global-local  model,  most  analytical  solutions  discussed  in 
this  section  are  obtained  by  using  various  engineering  methods  to  solve  the 
displacement-equilibrium  equations  under  certain  assumptions.  Thus,  these  can  be 
regarded  as  approximate  solutions  based  upon  elasticity  theory. 

2.2.1  Approximate  Elasticity  Solution 

Investigations  of  the  f  ree  edge  problem  ex  as  carried  out  by  Puppo  and  Hvensen  [2] 
using  a  composite  model  essentially  consisting  of  a  set  of  anisotropic  layers  separated 
by  isotropic  adhesive  layers.  It  was  assumed  that  the  isotropic  layers,  developed  only 
interlaminar  shear  stresses,  acting  as  an  adhesive  between  the  anisotropic  layers.  It  was 
reported  that  a  sharp  rise  of  the  interlaminar  shear  stress  could  be  observed  in  finite 
width  laminates.  However,  the  simplicity  of  these  elastic  formulations  precluded 
calculation  of  the  transverse  normal  stress,  and  the  problem  became  more  complicated 
w'hen  more  layers  were  involved. 

In  an  attempt  to  approximate  the  interlaminar  normal  stress,  a  simplified  formula 
was  developed  bv  Pagano  and  Pipes  [  1  ].  The  strategy  was  to  use  solutions  along  the 
longitudinal  mid  plane  ol  the  laminate  based  upon  classical  laminated  plate  theory,  one 
could  then  compute  the  force  and  moment  resultants  caused  by  the  interlaminar  stresses 
on  any  plane  z=constant  through  consideration  of  static  equilibrium.  The  maximum 
interlaminar  normal  stress  at  the  free-edge  could  then  be  expressed  in  terms  of  the 
transverse  stress  in  the  y-direction  calculated  from  the  laminated  plate  theory. 

Another  approximate  elasticity  solution  proposed  by  Pipes  and  Pagano  [3]  w-as  based 
upon  displacement-equilibrium  equations  for  an  anisotropic  elastic  medium.  Assuming 
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t he  transverse  stresses  in  the  y-,  z-  directions  to  \unish.  the  equations  were  written  in 
terms  of  the  single  variable  L  (axial  displacement  function).  This  yielded  components 
of  displacement,  strain  as  well  as  remaining  stress  fields  in  the  form  of 
sinusoidal-hyperbolic  series.  However,  violation  of  stress  equilibrium  in  the  transverse 
directions  as  w'ell  as  neglect  of  the  interlaminar  normal  stress  constituted  major 
drawbacks  of  this  scheme. 

2.2.2  Modified  Higher  Order  Theory 

1’aguno  [4]  derived  another  approximate  method  for  determination  of  distribution  <>t 
the  interlaminar  normal  stress  along  the  mid  plane  of  a  symmetric,  finite  width 
laminate.  The  approach  was  based  upon  a  modified  version  of  a  higher  order  theory 
proposed  by  Whitney  and  Sun  [28],  which  recognized  the  effect  of  shear  deformation 
through  the  inplane  rotations  as  w'ell  as  the  thickness  strain  implemented  in  the 
assumed  displacement  field.  However,  like  the  approximate  theories  discussed 
previously,  none  of  them  were  able  to  determine  the  complete  stress  field  near  the 
free-edge. 

2.2.3  Galerkin's  Method 

Due  to  the  fact  that  high  stress  gradients  occurring  near  the  free-edge  are  difficult 
to  estimate  by  numerical  approaches,  Wang  and  Dickson  [5]  applied  the  extended 
Galerkin’s  approach,  in  which  interlaminar  stresses  and  displacements  of  each  layer 
satisfying  geometrica1  boundary  conditions  were  represented  as  series  of  Legendre 
polynomials.  The  final  solution  was  reached  by  requiring  the  satisfaction  of  continuity 
conditions  at  each  interface  as  well  as  stress  boundary  conditions  at  exterior  planes. 
Due  to  tl  completeness  of  Legendre  polynomials,  convergence  of  solutions  could  be 
expected. 
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2.2.4  Perturbation  Technique 


In  an  effort  to  obtain  more  accurate  Iree-edge  stress  intensities,  a  perturbation 
technique  was  applied  by  Hsu  and  Herakovich  [6]  to  solve  the  three  coupled 
dimensionless  partial  differential  equations  based  upon  a  displacement  formulation  of 
the  elastic  problem.  It  showed  that  the  perturbation  solution  provided  a  smooth 
continuous  stress  distribution  in  the  vicinity  of  the  free-edge.  However,  this  solution 
had  the  limitation  that  the  shear  stress  distribution  was  a  function  of  both  laminate 
thickness-to-width  ratio  and  a  problem-dependent  parameter.  Although  the  latter  could 
be  chosen  such  that  the  maximum  values  of  shear  stress  field  did  not  exceed  elastic- 
limits,  the  accuracy  of  the  calculated  stresses  was  suspect. 

2.2.5  Boundary  Layer  Theory 

A  boundary  layer  theory  for  laminated  composites  in  plane  stress  was  developed 
by  Tang  and  Levy  [7]  from  the  three-dimensional  theory  of  anisotropic  elasticity.  By 
expanding  the  stresses,  displacements,  body  forces  and  surface  tractions  in  power  series 
of  the  half-thickness  of  a  lamina  in  the  equations  of  equilibrium,  compatibility  and 
boundary  conditions,  a  sequence  of  systems  of  equations  was  obtained.  The  complete 
solution  was  obtained  by  combining  solutions  of  the  interior  domain  based  on  the 
classical  lamination  theory  and  those  from  boundary  layer  and  matching  a  set  of 
appropriate  boundary  conditions.  This  formulation  indeed  provided  a  way  to  obtain 
analytical  solution  for  estimating  interlaminar  normal  as  well  as  shear  stress 


distribution. 


2.2.6  Reissner's  Variational  Principle 

In  order  to  have  displacement  as  well  as  stress  continuity,  a  mixed  formulation  is 
sometimes  used.  Unlike  the  elastic  approximations  discussed  previously,  Pagano  [8] 
developed  an  approximate  theory  for  a  general  composite  laminate  based  upon  an 
application  of  Reissner's  variational  principle.  In  this  theory,  the  inplane  stresses  are 
considered  linear  in  the  thickness  coordinate  while  the  transverse  stresses  derived  from 
equilibrium  consideration  are  cubic.  If  a  laminate  or  a  single  lamina  is  viewed  as  an 
assembly  of  N  sheets,  each  having  a  finite  thickness  and  required  to  satisfy  force  and 

moment  equilibrium,  the  analysis  led  to  a  set  of  23N  algebraic  and  ordinary 

differential  equations  which  had  to  be  solved  simultaneously.  Based  upon  the 

assumption  that  the  stress  field  is  independent  of  the  longitudinal  axis,  Pagano  [9] 
further  specialized  the  theory  to  the  free-edge  problem  by  reducing  the  stress  field 
determination  to  the  solution  of  a  one-dimensional  problem.  Despite  the  relative 
accuracy  of  this  theory  resulting  from  the  improvement  of  smoothness  for  both 
displacement  and  traction  fields  at  interfaces  between  adjacent  layers,  a  major  drawback 
was  that  its  application  was  limited  at  most  to  six  sublayers. 

2.2.7  Global-local  Model 

Pagano  [10]  introduced  a  global-local  model,  which  was  able  to  define  detailed 
response  functions  in  a  particular,  predetermined  region  of  interest  while  representing 
the  remainder  of  the  domain  by  effective  properties,  that  reduced  the  number  of 
variables  in  a  given  problem.  In  this  model,  for  the  global  region  of  the  laminate, 

potential  energy  has  been  utilized,  and  the  displacement  components  were  based  upon 

the  assumption  given  by  Whitney  and  Sun  [28].  The  Reissner  variational  principle 

described  in  [8],  however,  was  applied  for  the  local  region  in  which  a  thickness 
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distribution  of  stress  field  satisfying  equilibrium  equation  within  each  layer  was 
assumed.  A  variational  principle  was  then  used  to  derive  the  governing  equations  of 
equilibrium  for  the  whole  system.  It  was  reported  that  the  global-local  model  could 

effectively  solve  the  same  class  of  free-edge  stress  problem  as  described  in  [8]  and  had 

wider  range  of  applicability. 

2.3  Finite  Difference  Method 

2.3.1  Pesudo  Two-dimensional  Analysis 

Pipes  and  Pagano  [  1 1  ]  used  the  classical  theory  of  linear  elasticity  to  formulate  the 
problem  of  free-edge  delamination  of  a  strip  under  uniform  axial  strain.  Allowing  for 

material  symmetries  and  uniform  extension,  the  transverse  components  of  displacement 

were  assumed  to  be  independent  of  the  longitudinal  coordinate.  The  three  coupled 
elliptic  equations  for  the  displacement  functions  were  solved  using  a  finite  difference 
solution  technique  to  approximate  the  interlaminar  stresses.  Delamination  was  assumed 
to  be  primarily  due  to  the  high  shear  stress  near  the  free-edge  and  the  interlaminar 
stress  field  w'as  found  to  be  an  edge  effect  which  was  restricted  to  a  boundary  region 
approximately  equal  to  the  laminate  thickness. 

2.3.2  Three-dimensional  Analysis 

A  three-dimensional  finite  difference  analysis  was  carried  out  by  Altus,  Rotem  and 
Shmueli  [12]  to  examine  the  free-edge  stress  field.  The  displacement  equilibrium 
equation  was  solved  by  using  central  difference  method  while  for  displacement  or 
traction-free  boundary  conditions  as  well  as  interfacial  continuity  conditions,  either 
forward  or  backward  difference  scheme  was  applied.  Convergence  of  the  solution  was 
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expected  providing  a  reasonable  displacement  field  was  assumed  initially.  Although  a 
complete  stress  field  was  available  due  to  three-dimensional  characteristics,  an  iteration 
scheme  could  lie  a  serious  inconvenience. 


2.4  Finite  Element  Method 

In  order  to  more  effectively  evaluate  the  high  gradient  stress  field  at  the  free-edge 
of  laminated  composites,  the  popular  finite  element  method  has  been  applied  bv 
numerous  investigators. 

2.4.1  Displacement  Method 

Wang  and  Crossman  [13]  used  a  very  line,  constant  strain  triangular  element  grid 
to  model  the  laminate  boundary  region  through  a  cross-section.  The  functional 
dependence  of  the  assumed  displacement  field  was  of  the  same  type  as  in  Pipes  and 
Pagano's  analysis  [ll].  To  overcome  the  difficulty  of  computational  storage  and  time 
limitation,  the  solution  process  adopted  the  so-called  ’’sky-line"  matrix  storage  scheme. 
The  results  indicated  that  the  interlaminar  as  well  as  inplane  stress  singular  behavior 
was  highly  localized  in  angle-ply  laminated  composite.  A  simplified  method  for 
calculating  interlaminar  stress  was  proposed  [  1 4]  wherein  the  stresses  at  the  desired 
layer-interface  were  evaluated  by  substructuring  the  laminate  with  few'er  numltei  of 
effective  layers.  This  reduced  the  numbeT  of  laminar  interfaces  and  facilitated  linite 
element  calculation  within  fewer  elements. 

A  quasi-three-dimensional  finite  element  analysis  was  carried  out  by  Raju  and 
Crew  [15]  using  eight-noded  isoparametric  elements.  In  order  to  approximate  the  stress 
singularities,  polar  mesh  was  introduced  near  the  intersection  of  interface  and  free-edge, 
associated  with  a  so-called  log-linear  procedure  to  relate  the  steep  gradient  stress  with 
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the  radial  distance  from  the  singular  point  in  the  logarithmic  coordinate.  One  major 
drawback  of  this  scheme  is  that  the  power  of  singularity  has  to  be  determined  by 
solutions  calculated  from  finer  polar  mesh  near  the  interface  of  the  free-edge. 

Whitcomb,  Raju  and  Goree  [l  6]  further  pointed  out  that  the  disagreement  for  both 
magnitude  and  sign  of  the  interlaminar  normal  stress  distribution  among  various 
numerical  methods  could  be  attributed  to  the  unsymmetric  stress  tensor  at  the 
singularity.  In  their  approach  too,  the  problem  was  modeled  by  eight-noded 
isoparametric  elements.  It  was  concluded  that  finite  element  displacement  models  were 
capable  of  giving  accurate  stress  distributions  everywhere  except  in  the  region  within 
two  elements  of  a  stress  singularity. 

In  summary,  we  observe  that  in  the  conventional  displacement-based  finite  element 
formulation,  evaluation  of  shear  as  well  as  normal  stresses  required  expensive  mesh 
refinement  near  the  boundary  region  to  approximate  the  singular  stress  field.  Even 
then  the  actual  stress  distribution  along  the  free-edge  was  generally  not  sufficiently 
accurate. 

2.4.2  Stress  Method 

Rybicki  [17]  used  a  three-dimensional  equilibrium  finite  element  analysis  procedure, 
based  upon  minimization  of  complementary  energy,  to  solve  the  free-edge  stress 
problem.  Due  to  the  fact  that  the  assumed  stress  state  in  the  analysis  did  not  contain 
singular  term,  a  finite  rise  in  interlaminar  normal  and  shear  stresses  near  the  interface 
corner  was  observed  for  angle-ply  layup.  However,  this  method  involved  very  large 
matrices  and  was  computationally  expensive,  and  even  at  that  did  not  yield  a 
continuous  stress  field. 
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2.4.3  Hybrid  Assumed  Stress  Mode] 


In  Plan's  hybrid  model  [29],  stress  equilibrium  in  the  interior  of  the  elements  as 
well  as  displacement  continuity  along  interelement  boundaries  are  ensured,  but  the 
interelement  stress  continuity  is  satisfied  only  in  a  weighted  integral  sense.  Following 
Pian's  formulation,  Spilker  (18]  developed  a  special  hybrid  element  for  the  edge-stress 
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procedure  to  angle-pi v  laminates  because  in  these  the  complete  stress  field  has  to  be 
considered. 

A  special  formulation  of  a  singular  composite-edge  element  was  developed  by 
Wang  and  Yuan  [19]  based  on  the  Boundary  ’  -er  theory  [30]  and  the  variational 
principle  of  a  modified  hybrid  functional.  In  the  analysis,  the  singular  hybrid  element 
W'as  used  in  conjunction  with  displacement-based  eight-noded  isoparametric  elements,  and 
it  was  reported  to  give  satisfactory  stress  distribution  near  the  free-edge.  This  method 
is  excellent  for  determining  possible  growth  of  delamination  but  would  be  awkward  to 
use  to  predict  occurrence  of  delamination  in  an  intact  specimen.  This  is  because 
sometimes,  it  is  hard  to  find  the  nlace  in  which  stress  singularity  may  occur. 


2.5  Summary  and  Research  Motivation 

The  analytical  and  numerical  solutions  discussed  above  for  the  free-edge  stress 
problem  are  summarized  in  Table  ( 1 ).  Some  conclusions  can  be  made  at  this  point. 
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1.  Generally  speaking,  an  analytical  solution  lor  the  complete  stress  field  is 
extremely  difficult.  The  solution  procedure  for  the  case  of  a  multi-layer 
system  is  not  currently  available. 

2.  Use  of  a  finite  difference  technique  suffers  from  geometric  limitations. 

Calculation  of  stresses  at  the  interfaces  or  laminate  boundaries  needs  to  apply 
additional  techniques,  such  as  iteration  scheme.  Even  then  the  solution 

generally  lacks  credibility. 

3.  Conventional  displacement  based  finite  element  methods  are  incapable  of 

predicting  accurate  stress  fields  particularly  along  element  boundaries.  Stress 
equilibrium  approach  is  apparently  impractical.  Use  of  hybrid  element  does 
improve  stress  calculation  but  is  applicable  only  to  some  sjtecial  cases. 
Application  of  singular  element  near  the  free-edge  boundary  apparently  makes 
the  analysis  too  subjective.  In  order  to  have  reliable  predictions  of 
displacement  and  stress  fields,  it  is  necessary  that  the  free-edge  stress  model  be 
able  to  approximate  the  real  life  situation  as  closely  as  possible.  In  other 
words,  the  displacement  and  stress  continuity  conditions  along  with 

traction-free  boundary  condition  have  to  be  exactly  satisfied.  Considering  a'so 
the  generality  and  effectiveness  of  the  analysis,  the  displacement-based  finite 
element  approach  with  higher  order  interpolation  function  could  conceivably  be 
superior  to  the  other  approximate  theories. 
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Table  1:  Comparison  of  Various  Methods  for  Sohing  the  ITT)  Problems 


Ref 

Author 

method  of  analysis 

2 

Puppo  &  Evensen 

Elastic  approximation 

<7  O'  7  r 

>  y  \Z  xy 

3 

Pipes  &  Pagano 

Approximate  elastic  solution 

<J  7  7  7 

X  X\  \7  V/ 

4 

Pagano 

Modified  higher  order  theory 

<T 

7 

5 

Viang  &  Dickson 

Extended  Galerkin’s  approach 

o  t 

i  V7 

6 

Hsu  &  Herakovich 

Perturbation  technique 

(77  T 

Z  \7  \7 

7 

Tang  ic  Levy 

Boundary  Layer  theory 

<J  <J  <J  T  7  7 

\  y  z  y?  \?  x\ 

8 

Pagano 

Reissner's  variational 

principle- -mixed  method 

o'  a  a  r  7  7 

x  v  2  yz  \z  xy 

10 

Pagano  &  Soni 

Global-local  model 

a  (T  a  7  7  7 

s  y  z  yz  \z  xy 

11 

Pipes  &  Pagano 

Finite  difference  method 

(X  O'  (7  7  7  7 

x  y  z  yz  \7  xy 

13 

Wang  h  Crossman 

Finite  element  method: 
constant  strain  triangle 

(7  a  (7  7  7  7 

x  y  z  yz  \7  xy 

16 

Whitcomb  et  al . 

Finite  element  method: 
8-noded  isoparametric 

element 

a  a  a  t  7  7 

\  v  2  w  \v  x\ 

17 

Rybicki 

Finite  element  method: 
equilibrium  stress  approach 

(7  (7  (7  7  7  7 

x  y  z  yz  \z  \> 

18 

Spilker 

Finite  element  method: 

hybrid  assumed  stress  model 

o'  a  7 

y  z  yz 

19 

Wang  &  Yuan 

Finite  element  method: 
Singular  hybrid  element 

O’  O'  O’  7  7  7  ! 

x  y  z  yz  \z  x\ 

1ft 


SECTION  Ill 


VARIATIONAL  FORMULATION  AND  FINITE 
ELEMENT  APPROXIMATION  IN  LINEAR  ELASTICITY 

3.1  Introduction 

In  this  section,  a  variational  formulation  of  three-dimensional  elasticity  is  described 
and  its  use  as  the  basis  of  a  finite  element  approximation  is  discussed.  The  treatment 
essentially  follows  that  in  reference  [31].  Variational  formulation  has  been  used  as  the 
basis  lor  direct  methods  of  obtaining  approximate  solutions  to  boundary  value  and 
initial  boundary  value  problems.  Traditionally,  the  approximation  space  is  generated  by 
complete  orthonormal  sets  consisting  of  eigenfunctions  of  self-adjoint  operators.  The 
functions  which  are  used  to  approximate  the  field  variables  are  required  to  satisfy 
certain  continuity  requirements  over  the  whole  domain.  The  finite  element  method, 
however,  offers  an  alternative  route  for  generating  the  sequence  of  finite  dimensional 
approximation  spaces.  The  region  under  consideration  is  subdicided  into  a  finite 
number  of  discrete  elements,  and  the  field  variables  are  represented  by  functions  which 
follow  the  same  continuity  condition  only  piecewise  within  each  element.  Some 
significant  differences  between  the  finite  element  method  and  the  traditional  direct 
methods  include  [31] 

1-  The  base  functions  have  local  support  and  are  nonorthogonal. 

2.  The  sequence  of  approximation  spaces  is  ordered  by  refinement 

3.  The  local  support  functions  may  have  only  limited  smoothness 
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The  support  of  eucli  base  function  is  confined  to  the  neighborhood  of  a  nodal 
point  and  extends  over  the  elements  of  the  finite  element  approximation  sharing  that 
point.  Across  interelement  boundaries  within  the  support  and  at  the  boundary  of  the 
support,  the  function  may  have  only  limited  smoothness.  In  a  sequence  of  refinements, 
additional  nodal  points,  elements  and  base  functions  are  introduced.  The  base  functions 
associated  with  each  nodal  point  change  with  refinement,  including  a  monotonic 
decrease  m  support.  Additional  discontinuities  might  be  introduced  at  each  refinement. 
Variational  formulations  and  solution  procedure  for  direct  methods  based  on  the  finite 
element  approach  must  allow  for  these  peculiarities  of  finite  element  approximation 
spaces. 


3.2  Boundary  Value  Problem 

Consider  an  open  connected  region  R  in  an  euclidean  space.  QR  is  the  boundary 
of  R  and  R  its  closure.  A  typical  boundary  value  problem  on  R  is  defined  by  the  set 
of  equations 


Au  =  f  on  R 

Cu  =  g  on  $R 

where  A  is  the  field  operator  and  C  is  the  boundary  operator  such  that 

A:»irV„ 

C:,)*rv* 


(1) 

(2) 

(3.i 

(4) 


V r,  \'aK  are  linear  vector  spaces  whose  elements  are  defined  on  the  regions  indicated 
by  the  subscripts.  DR,  DaR  are  dense  subsets  in  VR,  VdR,  and  denote  the  domains  of  A, 
C  respectively.  ffR  is  the  extension  of  DR  i.e.  any  element  u€DR  has  a  unique 
extension  in  D;iR  and  every  element  in  l)oR  is  the  extension  of  an  element  (not 


necessarily  unique)  in  DR.  For  given  f€VR,  g€V9R,  the  boundary  value  problem 
consists  of  determining  u€DR  along  with  its  extension  in  DSR  such  that  (l)  and  (2)  are 
satisfied. 

3.3  A  Variational  Principle 

Let  the  linear  operator  A  be  self-adjoint,  i.e.  there  exists  a  nondegenerate,  linear 
Gateaux  differentiable,  bilinear  mapping  BR:DRxVR-*S,  where  S  is  a  linear  vector  space, 
such  that 

BR(u,Av)  =  BR(v,Au)  +  CaR(v,u)  u,v6DRDVR  (5) 

Here,  CflR(v,u)  are  quantities  associated  with  the  boundary  QR.  Magri  [32]  has  shown 
that  such  a  bilinear  mapping  can  be  constructed  for  every  linear  operator  A.  If  the 
boundary  operator  C  is  consistent  [33]  with  the  field  operator  A,  i.e.,  there  exists  a 
nondegenerate,  linear  Gateaux  differentiable,  bilinear  mapping  Bg^D^xYg,,-^,  such  that 

C5R(v,u)  =  BdR(v,Cu)  -  B^Cu.Cv)  (6) 

then,  the  linear  Gateaux  differential  of 

ft(u)  =  BR(u,Au-2f)  +  B8R(u,Cu-2g)  (7) 

vanishes  if  and  only  if  (l),  (2)  are  satisfied.  Sandhu  and  Salaam  [33]  further  pointed 
out  that  even  if  the  boundary  condition  is  homogeneous,  i.e.  g  =  0,  the  quantity 
BdR(u,Cu)  in  (7)  must  be  included  if  the  variational  principle  is  to  hold  for  the  path 
of  Gateaux  differentiation  not  satisfying  homogeneous  boundary  conditions.  This  is 
important  for  approximation  in  finite  element  spaces  where  the  variation  is  introduced 
as  change  in  the  nodal  point  value  of  the  field  variable  and,  consequently,  the  path  of 
variation  may  not  satisfy  the  boundary  condition  or  internal  smoothness. 
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3.4  Variational  Principle  for  Finite  Element  Approximation 


In  the  finite  element  method,  the  region  R  is  approximated  by  a  set  of  elements 
]Re;  e=l,2 . m}  such  that 

ReflRf=0  if  e*f  (8) 

m 

lim  [JR  =  R  (9) 

ID  —  oo  ^  * 

l*=l 

The  field  variables  are  approximated  hv  I  unctions  which  may  not  be  sufficiently 
sm<x>th.  However,  oxer  each  element,  adequate  smoothness  is  assured.  If  R  represents 
the  interior  of  the  e-th  element  ;md  QR  its  ktundarv.  we  have  [34] 

B  (u,Av)  =  B  (v,Au)  +  C  ,  (v,u)  (10) 

K  f*K 

and 

CaR  (v,u)  =  BaR  (v.Cu)  -  BaR  (u,Cv)  (11) 

Further  define 

in 

n(u)=Z[BR(u.Au-2f)  +  B  n  (u,Cu-2g)]  +  B  ,(u,(Cu)')  (12) 

e=i  e  e  0K* 

where  QR[  represents  interior  boundaries  of  elements  and  a  prime  denotes  a  jump.  The 
Gateaux  differential 


A  Mu)  =  2  T[B.,  ( v,Au-f)  +  R  ...  (v.('u-g)]  +  2  B  (v,(Cu)') 

*  Hh  I 


(13) 


e=l 


vanishes  if  and  only  if  (1),  (2)  are  satisfied  over  each  element  and  (Cu)1  vanishes,  i.e. 
Cu  is  continuous  across  interelement  boundaries.  If  there  are  actual  discontinuities  in 
the  interior  of  R,  let 

(Cu)'=g'  over  QR1  (14) 
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where  g'  is  specified  over  QR1.  Then,  if  the  union  of  intersection  of  element 
boundaries  covers  R1,  the  functional  in  (12)  may  be  redefined  as 

m 

ft(u)  =  £[Br  (u,Au-2f)  +  BdRneR  (u,Cu-2g)]  +  B  (u.(Cu)'-2g')  (15) 

e=  1  *  aR* 

3.5  Linear  Operator  with  Adjoint  Splitting 

Many  physical  problems  can  be  written  in  the  form  of 

Au  =  l  u  +  T  l:Tu  =  f  on  R  (16) 

where 

F;  WR-*V’R  (17) 

T :  WR-*XR  (18) 

E:Xr-Yr  (19) 

f:YR-VR  (20) 

T*  is  the  adjoint  of  T,  i.e.  BR,  BR  such  that 

Br(u,Tv)  =  Br(v,T  u)  +  CaR(v,u)  (21) 

Here  Br:WrxXr-»S  and  BR:WRxVR-*S.  S  is  a  linear  vector  space  and  BR,  BR  are 
continuous  non-degenerate  bilinear  mappings.  !■,  1:  are  symmetric,  i.e. 

BR(u,Fv)  =  Br(v,F.u)  (22) 

Er(u,Fv)  =  B^v^u)  (23) 

Introducing  €,  a  through  the  equations 

Tu  —  e  =  0  on  R  (24) 

Ee  -  cr  =  0  on  R  (25) 

(16)  can  be  written  as 
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Fu  +  T  ar  =  f  on  R 

Combining  (24)  through  (26),  these  constitute  the  coupled  system 


(26) 


F  0 

T 

u 

f 

0  E 

-1 

£ 

= 

0 

T  -1 

0 

a 

0 

on  R 


(27) 


If  the  inverse  of  I:  exists,  let  G=I:'.  Then,  combining  (24),  (25) 

Tu— (icr=0  on  R  (28) 

(26)  and  (28)  are  the  coupled  system 


1  T 

u| 

If 

T  -Cl 

0*1 

1" 

(29)  is  referred  to  as  the  complementary  form. 


For  an  operator  with  adjoint  splitting,  let  the  boundary  conditions  on  u,  cr  be 


Ciu  =  8. 

on  S,  0R 

(30) 

C2a:=&2 

on  S2  QR 

(31) 

The  discontinuity 

conditions  are 

(C,u)'  =  g” 

.  on  S', 

(32) 

(C2cr)'=g'2  on  S‘2  (33) 

where  S,  and  S'?  are  interior  surfaces  imbedded  in  the  intersection  of  finite  element 


boundaries.  C,,  C2  consistent  with  T,  T*  implies  the  existence  of  bilinear  mapping 
Bsf  such 

Bh  (u,Tv)  =  B  (v.Tu)  +  B  (v,C,u) - B  (u.C.v)  (34) 

Ke  Re  S*  2  S'  1 

where  S2,  S*  are  complementary  subsets  of  boundary  0Re  of  element  e. 
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The  function  governing  variational  formulation  of  (27)  for  finite  element 
approximation  is 

m 

n(u,e,tr)  =  £[Br  (u,Fu+T  cr— 2f)  +  BR  (e,Ee-tr)  +  BR  (cr.Tu-e)  +  B8R  ns  (o-,C1u-2g])] 

e=  1 

+  fBaRenS2(u,C2cr-2g2)]  +  BjCtr  ,(C,u)  — 2g” , )  +  B'2(u,(C2o-)'-2g'2)  (35) 

It  is  important  to  note  that  even  if  there  are  no  interior  discontinuities  in  the 
physical  problem  or  the  specified  boundary  conditions  are  homogeneous,  the  boundary 
and  the  discontinuity  terms  must  be  included  to  accommodate  the  nature  of  the  finite 
element  approximation  space  [32]. 


3.6  Principle  of  Minimum  Potential  Energy 


The  field  equations  for  isothermal  quasi-static  deformation  of  anisotropic,  linear 
elastic  solids,  assuming  no  initial  stresses  and  strains,  are: 
a)  Equilibrium  of  stresses 

cr  +f  =0  on  R  (36) 

>>>  j 


b)  Kinematics 

For  small  deformation,  the  strain-displacement  relationship  is 


u.  i  =  € 

(i.j)  ij 


on  R 

c)  Constitutive  relations 


cr  =  E  .  . 

ij  ijkl  kl 


(37) 


(38) 


on  an  open  bounded  connected  set  R  contained  in  the  three-dimensional  Euclidean  space 
E.  Here  u,,  f,,  €kl,  cr^,  E1Jkl  are,  respectively,  the  components  of  the  displacement  vector, 
the  body  force  vector,  the  infinitesimal  strain  tensor,  the  symmetric  Cauchy  stress 
tensor  and  the  isothermal  elasticity  tensor.  The  range  of  indices  is  1,  2,  3  and 
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summation  on  repeated  indices  is  implied.  A  subscript  following  a  comma  denotes 
partial  differentiation  with  respect  to  the  coordinate,  in  the  reference  frame,  defined  by 
the  subscript.  Parantheses  around  subscripts  denote  the  symmetric  part  of  the  quantity. 

Let  the  functions  ui(  €ij;  cr^  satisfy  the  continuity  and  differentiability  properties 
required  in  the  equations  of  elasticity  over  every  subregion  Rr  Then,  admitting 


(ui,€ij,cri.)  as  the  15-tuple  of  dependent  variables,  components  of  vectors  and  tensors 
being  regarded  as  ordered  subsets  in  an  n-tuple,  ( 36)-( 38)  can  be  written  as  [33] 


0 
0 

1(8  A+8  J_)  -i 

2  k'  61  11 0k 


0  — (S„-l+8  L-l) 

2  'k  0j  Jk  0i 


'ijkl 


-1 

0 


f 

k 

ek, 

= 

0 

cr 

0 

'j 

on  R 


(39) 


Consistent  boundary  conditions  for  the  problem  are 


— n  u.  =  — n  Q 

j  < 

on  S, 

(40) 

a..n.  =  l. 

on  S2 

(41) 

•j  '  j 

where  the  nj  are 

components  of  a  unit  normal 

to  the  boundary  S,  and  the  jump 

conditions  are 


(tr  n.)'  =  g'  .  on  S‘ 

lj  i  ®  2j  2 

— (nu)'  =  —  g'.  on  S1, 

j  I  ®  Ijl  i 

Setting  up  the  problem  in  inner  product  space,  i.e. 


(42) 

(43) 


uvdR,  and  defining 


br(u,v)  =  XrBR(u’v)R,  <44> 

e=  1 

The  basic  functional  corresponding  to  (35),  allowing  for  relaxed  continuity,  is  [33] 
ft  ,(u  ,€.  „cr  .)  =  f  (e  £  ...€—  ucr  — 2e  <r  -2u  f  +u  rr)dR 

I'  i’  kl’  i/  J  'j  ijkl  kl  1  1J.J  ij  1J  11  IJ 
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(45) 


+  / 

—  r  cr  .((n.u.)'— 2g' .  )dS 
J  >j  j  i  6  »  j* 
s 

In  using  (45)  as  the  basis  for  finite  element  approximations,  it  is  not  necessary  for  the 
interpolants  to  satisfy  any  boundary  conditions  or  interelement  continuity.  For  no 
jump  discontinuities,  g'tJ1  and  g'2j  vanish.  However,  it  is  important  to  retain  the  terms 
containing  (crijni),1  (mu^X  in  the  formulation. 

Using  symmetry  property  of  the  operator  matrix,  i.e. 


u  (or  n  —  2f)dS—  f  cr  n  (u  —  2u  )dS+  f  u  ((cr  n  )'— 2g,  )dS 
j  *j  >  j  J  s  ‘J  J  1  1  J  ,  J  IJ  1  2J 


/cr  u  dR  =  —  /  cr  u  dR  +  /  an  u  dS  +  /  cr  (n.u  )'dS 
R  ‘>J  '  J  R  ,J  ‘,J  J  S  ,J  J  '  J  ci  U  J  ‘ 

+  I*  uj(a..nj)'dS  (46) 

J  s‘ 

to  eliminate  the  term  containing  cr,^  from  (45),  the  functional  can  be  written  as 


ft  ,(u  ,e.  „cr  )  =  f  €.  E. .  €,  dR  +  2  /  cr  (u  — e  )dR  — 2  f  uf  dR 

2  kl’  «/  J  R  >J  .jkl  kl  J  R  U  ‘-J  ‘J  Jr  " 

—  2  f  ut  dS  — 2  f  cr  n  (u  — Q  )dS  —  2  f*  (r  (nu  )' 

J  s''  J  s'3  3  '  '  J  c  ,J  J  1 


(47) 


ft2  is  the  modified  variational  principle  with  three  independent  fields  proposed  by 
Prager  [35].  If  u,,  el}  are  restricted  to  satisfy  the  last  of  (39),  the  strain-displacement 
relations,  Prager's  modified  principle  of  total  energy  theory  is  obtained 

ft  =  f  €  E  6.  ,dR  —  2  f  u.f.dR  —  2  f  u.t  dS  —  2  f  cr  n  (u  — fi  )dS 
3  J  R  'J  'Jkl  kl  Jr"  J  s'  '  j  s'33'  ' 


—  2  f  cr  (n  u  )'dS 

J  ,J  J  1 


(48) 
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(48)  was  also  proposed  bv  Pian  and  'long  [29]  and  is  the  basis  of  their  hybrid  method 
with  assumed  displacement  field.  If  the  displacement  field,  u,,  is  further  restricted  to 
satisfy  the  displacement  boundary  condition  (40)  on  S,,  reduces  to 


n  =  f  €  E  €  dR  —  2  f  ufdR-2  f  \1tdS-2f  cr  (n  u  )'dS 

4  J  R  ‘J  ‘Jkl  kl  JR"  J  S'  '  J  c.  U  J  ' 


(49) 


If  the  finite  element  interpolations  are  chosen  to  identically  satisfy  displacement 
continuity  across  S',,  the  last  term  in  vanished  and  (49)  becomes 

ft  =  f  e  F,  e,dR-2  f  ufdR-2  f  u  f  dS  (50) 

Jr  Jr  J  s, 

The  vanishing  of  the  variation  of  Qs  with  respect  to  the  displacement  components  u, 
implies  the  satisfaction  of  equilibrium  equations  (36).  This  functional  corresponds  to 
the  classical  principle  of  minimum  potential  energy  which  is  customarily  used  as  the 
basis  of  the  finite  element  displacement  formulation  of  the  elastostatics  problems. 


3.7  Assumed  Displacement  Finite  Element  Formulation 

For  the  boundary  value  problem  stated  in  (l)  and  (2),  the  solutions  u  to  the 
forcing  functions  f  in  general  belong  to  L2,  the  space  of  square  integrable  function.  L2 
is  a  separable  Hilbert  space.  However,  u  may  be  contained  in  a  subset  D  of  L2  such 
that  A,  the  linear  operator,  is  defined  on  D.  We  assume  that  D  is  dense  in  L2.  If 

the  set  of  functions  {<£k,  k=l,2, oo}  is  a  basis  in  D,  then  any  function  u€L2  can  be 

expressed  as  an  infinite  sum: 

OO 

U  =  2>A  (51) 

k  =  1 
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A  scheme  to  generate  approximate  solutions  is  to  use  a  finite  set  of  terms  in  the 
infinite  sum  above.  Thus,  as  an  approximation 


u  =  £ak<*>k  (52) 

k  =  1 

The  approximation  process  consists  of  an  appropriate  choice  of  n,  <f>k  and  the  coefficient 
ak.  Several  alternative  procedures  are  available.  The  finite  element  method  is  a 
special  process  of  selection  of  finite  subset  of  the  basis  {$J.  The  coefficients  ak  are 
generally  evaluated  by  requiring  the  approximate  solution  to  satisfy  a  variational 
principle. 

The  finite  element  idealization  essentially  partitions  the  spatial  region  R  into  a 
finite  number  of  nontrival  discrete  elements  or  subregions.  The  geometry  of  the 
elements  is  defined  by  a  set  of  points  in  space  called  the  nodal  points  of  the  system. 
Over  an  element  m,  let  an  approximation  to  u  be  u"  such  that 


n 


U 


m 

n 


z<x 


k=l 


(53) 


or  in  matrix  form,  dropping  the  subscript  n, 

um  =  {<nT{am}  (54) 

where  {^m}1  is  a  row  vector  consisting  of  0k  as  its  elements  and  {ara}  is  a  column 
vector  of  coefficients  ak.  Evaluating  the  function,  and  its  derivatives  up  to  a  certain 
order  at  nodal  points,  yields 

{uim}  =  [^im]T{am}  (55) 

where  {uj"}  is  the  vector  of  nodal  point  values  of  the  function  and  its  derivatives  up 
to  the  order  selected,  and  is  the  matrix  of  base  functions  evaluated  at  each  nodal 
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point.  The  rows  and  columns  of  are  linearly  independent.  If  square,  the  matrix 

is  invertible.  Hence,  we  can  write 

{am}  =  ([<£im]T)  1  {u^}  =  [A]  1  {u"'}  (56) 

where  A=[$"Jr 

Substituting  (56)  into  (54) 

um=[rf[Ar'{u;"}  =  [0m]T{u™}  (57) 

where  [0ra]  can  now'  be  regarded  as  a  set  of  interpolation  functions  relating  nodal  point 
values  of  a  function  and  its  derivatives  up  to  a  preselected  order,  to  the  value  at  an 
arbitrary  point  within  the  element  m  [36]. 

In  applying  the  potential  energy  functional  shown  in  (50),  is  assumed  to  satisfy 
the  strain-displacement  relationship,  and  the  displacement  field  should  satisfy  the 
prescribed  displacement  boundary  condition  on  S,.  Vanishing  of  the  variation  would 
imply  satisfaction  of  the  equilibrium  equations.  As  stated  previously,  in  the  finite 
element  method,  the  displacement  u  is  approximated  by  interpolation  functions  and 
generalized  displacements  at  a  finite  number  of  nodal  points  of  each  element.  The 
interpolation  function  must  be  chosen  in  such  a  way  that  when  the  nodal  point 
displacements  for  two  adjacent  elements  are  compatible,  the  displacements  along  the 
common  boundary  are  compatible.  Meanwhile,  the  interpolation  function  must  also 
satisfy  the  requirement  that  the  first  derivatives  of  the  displacement  field  exist. 

Based  on  (57),  the  assumed  displacement  over  an  element  can  be  rewritten  in 
matrix  form  as 


{un5}=[0qml^]T 


(58) 
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where  qm  is  the  column  matrix  of  generalized  displacements  at  boundary  nodes  of 
element  m  which  determines  the  inter-element  continuity,  and  rm  is  the  column  matrix 


of  generalized  displacements  at  nodal  points  either  at  the  boundary  or  at  the  interior 
of  element  m  but  which  do  not  affect  the  interelement  continuity  requirements.  The 
corresponding  strain  distribution  is 

ml 


{en,}  =  (OC)T 


(59) 


where  <5™  and  0™  are  obtained  by  differentiating  0"  and  0"'  with  respect  to  the 
spatial  coordinates.  Substituting  (59)  into  (50)  and  expressing  in  the  finite  element 
discretized  form,  we  have  [37] 


M  „  I 


[Cq'OEmH0;qi  oTl 


1 

dRm~2  J  *  [0'n!0;'][0fm]T{f}dRo 


-»/ 


s,na» 


(60) 


where 


[Em]=  matrix  of  elastic  constants  for  element  m 
[0™]=  matrix  of  interpolation  functions  for  the 
body  forces  in  element  m 
[0™3=  matrix  of  interpolation  functions  for  the 

prescribed  tractions  on  the  surface  Sm  of  element  m 

The  summation  sign  in  (60)  implies  the  direct  stiffness  assembly  procedure,  and  the 


vector 


is  the  vector  of  global  displacements. 


fis  can  also  be  written  as 


O 5  =  Z«qm>K] {qm)  +  2  {rm}T[K*] {qm}  +  (rm}T[K^] {rm} 

m  =  l 


-2{Fqm}T{qm}-2{F^V}) 


(61) 
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where 


K  =  /  [0m][F£m][0m]TdR 

qq  I  LYeqJL  1  L^eqJ  m 

m 

(62) 

K  =  f  l4>m][Em)[4>m?dR 

rq  J  lTer  1  A  m 

m 

(63) 

K  =  f  [4>m][En'][<Am]TdR 

rr  I  ”eq  1  JLTerJ  m 

J  R 

m 

(64) 

Fq=  J  [0'qn][0jn]T{f,n}dRiii+  J 

m 

r  i</>"‘3T{tm}dsm 

1  sns... 

(65) 

f  =  f  [0;"][0;"]T{f"'}dRn]  + 

•i  R  « 

m 

f  [0;"][<f{tm}dsm 

'  s2nsm 

(66) 

The  displacements  {rm}  in  element  m  are  independent  of  displacements,  {r1},  for  i^m. 
The  stationarity  condition  with  respect  to  their  variations  yields 

[K™]  {qm}  +  [K™]  {r"1}  —  {F™}  =  0  (67) 

Solving  (67)  for  {rm}  yields 

{rm}  =  [K™r1({F™}-[K™]{qm})  (68) 

Substituting  (68)  into  (61)  yields 

M 

=  £({cT}T[Km:i  {qm}  -  {Fm}T{qm}  +  Cm)  (69) 

m  =  l 

where  [Km]  and  {Fm}  are,  respectively,  the  element  stiffness  matrix  and  the  equivalent 
nodal  forces  defined  by 

[Km]  =  [Kqmq]  -  [K^/lK^r  ‘[K™]  (70) 

{Fm}  =  {Fq  }  -  ‘{F™}  (71) 

Cm  =  — {F™}  [K  ™F*{F™}  =  constant  (72) 
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fts  in  (69)  is  given  in  terms  of  the  generalized  displacements  {q}  which  are  not 
independent  for  different  elements.  Using  global  coordinates,  (69)  can  be  written  as 

n  5  =  {q}T[K]  {q}  -  2  {q}T{F}  +  C  m  (73) 

Taking  the  variation  of  this  discretized  form  of  the  functional  yields  the  system  of 
algebraic  equations 

[K]{q}  =  {F}  (74) 

which  can  be  solved  for  the  unknown  nodal  displacement  (q).  The  matrix  [K]  is 
positive  definite,  symmetric  and  banded.  The  process  of  eliminating  the  generalized 
coordinates,  {r},  from  each  element  is  called  the  static-condensation  process  [38].  The 
introduction  of  these  terms,  which  do  not  interfere  with  interelement  compatibility, 
results  in  an  improvement  in  the  satisfaction  of  the  equilibrium  equations  within  each 
element.  However,  the  satisfaction  of  the  equilibrium  equations  along  the  interelement 
boundary  is  still  governed  by  the  degrees  of  compatibility  supplied  by  the  interpolation 
functions  for  the  generalized  displacements,  {q}.  The  solution  obtained  represents  an 
underestimate  of  the  true  solution  in  the  sense  of  energy  [39]. 
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SECTION  IV 


CONTINUOUS  STRAIN  FINITE  ELEMENT 
INTERPOLATION 


4.1  Introduction 

In  the  finite  element  method,  the  displacement  field  is  approximated  b\ 
interpolation  functions  and  generalized  displacements  at  a  finte  number  of  nodal  points 
which  also  define  the  geometry  of  the  elements.  To  ensure  continuous  strain  across 
interelement  boundaries,  it  is  sufficient  that  the  interpolation  functions  be  such  that 
the  displacement  components  as  well  as  their  first  derivatives  along  the  common 
boundary  are  continuous. 

Tocher  and  Hartz  [40]  pointed  out  that  for  plate  bending  analysis,  continuity  of 
slopes  of  the  plate  displacement  surface  is  necessary.  The  compatible  cubic  interpolation 
functions  developed  by  Tocher  [40]  and  by  Clough  and  Felippa  [4 1  ],  among  others, 
satisfy  this  requirement.  For  plate  bending,  the  generalized  displacements  used  were  w, 
the  transverse  displacement  of  the  plate  and  its  derivatives  ws,  wy.  Mere,  the 
subscripts  x  or  y  denote  partial  differentiation  with  respect  to  the  independent 
variables  x,  y.  Applying  the  same  displacement  interpolation  scheme  to  the  plane 
elasticity  problems  [40],  the  corresponding  generalized  displacements  were 
u,  ux,  uy,  v,  vx,  vy,  the  in-plane  displacements  and  their  first  derivatives  at  each  node. 
Thus,  continuity  of  strain  between  adjacent  elements  was  ensured. 
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Instead  of  using  a  local  ('artesian  coordinate  system  in  the  derivation  of  the  cubic 
polynomial  with  nine  coefficients  for  each  displacement  component  in  Tocher  and 
Hartz's  work,  triangular  coordinates  were  used  in  the  present  fully  compatible 
quadrilateral  element,  following  Felippa's  [42]  work  on  plate  bending  analysis.  This 
simplified  the  generation  of  various  matrix  relationships  for  the  constituent  triangular 
elements.  Tocher's  [40]  element  used  incomplete  cubic  polynomials.  Felippa's  elements 
were  based  on  complete  cubic  polynomials  and  were  therefore  selected  for  application  to 
the  Iree-edge  problems.  These  elements  include  ’I. teller's  formulation  as  specialization. 
Cubic  expansion  of  the  9 degree  of ■  freedom  conforming  triangular  element  (LOCI'-*))  for 
both  in  plane  displacement  components  as  used  by  Tocher  [40]  was  extended  to 
quadrilateral  element  designated  (.1-15.  Quadrilateral  elements,  Q-19  and  Q-23, 
assembled  from  LCCT-11  and  LCCT-12  triangular  elements  introduced  by  Felippa  [42], 
were  also  redeveloped  for  the  plane  elasticity  problems.  The  continuous  strain  elements 
were  used  to  analyze  a  pseudo  two-dimensional  free-edge  stress  problem  similar  to  that 
of  Pipes  and  Pagano  [ll]  for  composite  laminate  coupons  under  uniform  extension. 

4.2  Interpolation  Functions  of  Continuous  Strain  Elements 

In  the  following.  Felippa’s  [41.42]  approach  for  deriving  the  plate  bending 
interpolation  functions  is  summarized.  We  use  the  same  element  name  as  Felippa's  and 
start  with  u  instead  of  w  for  the  plane  elasticity  problems.  Similar  derivations 
applied  to  the  displacement  component  v. 


4.2.1  LC'CT-12  Element 


A  complete  cubic  polynomial  in  two  variables  is  defined  by  ten  independent 
coel  1  icients.  1  he  values  ol  u,  the  x -direction  displacement  of  the  plane  stress  body  and 
its  derivatives  ux,  uv  at  the  three  vertices  of  a  triangle  yields  nine  independent 
quantities.  To  ensure  continuity  of  derivative  ur  across  element  boundaries,  it  is 
necessary  that  un  Ite  known  at  some  points  other  than  the  vertices  along  each  of  the 
three  edges.  It  is  convenient  to  introduce  mid  side  nodes  on  each  of  the  three  edges. 
1  his  element  then  has  twelve  independent  quantities  against  the  minimum  of  ten 
needed  to  completely  define  a  cubic  polynomial. 

In  order  to  use  a  cubic  polynomial  with  continuous  first  derivatives  in  the  interior 
as  well  as  on  the  element  boundaries.  Felippa  proposed  that  the  element  be  made  up  of 
three  subtriangles  as  illustrated  in  Figure  (2).  Kach  subtriangle  has  three  vertices  and 
one  mid-side  node  to  supply  the  ten  independent  quantities  for  defining  the  cubic 
polynomial  interpolation  in  its  interior.  The  point  C)  could  be  any  interior  point. 
However,  for  simplicity  of  formulation,  the  centroid  is  generally  used  [41 J. 

The  nodal  displacement  degrees  of  freedom  to  be  considered  in  the  stiffness  matrix 
of  the  complete  element  (Figure  2)  include  the  values  of  the  in-plane  displacement 
components,  u,  v.  along  with  their  first  derivatives  u  ,  u,  ,  v  ,  \  (i  =  1,2,3)  about  the 

x  and  y  axes  at  each  corner  as  well  as  the  normal  slopes  at  the  three  mid  side  nodes 
about  axes  perpendicular  to  these  sides  respectivel v,  viz.  u  .  u  ,  u  ,  and  v  ,  v  ,  v  . 

After  forming  the  expression  of  the  cubic  displacement  patterns  in  the  three 
subelements,  because  of  the  common  displacements  imposed  at  the  nodes,  the  in  plane 
displacements  of  two  adjacent  subelements  are  identical  along  their  juncture  line.  To 
establish  continuity  of  ur  along  the  edges  of  the  subelements,  it  is  sufficient  that  ur 


evaluated  at  points  7,  8,  V,  mid  points  ol  these  edges.  !  iom  adjacent  subtriangles  be  the 
same.  I'liese  three  conditions  were  used  to  evaluate  the  in-plane  displacement  u,  and 
its  derivatives  ux,  u  at  the  interior  point  ().  With  the  interior  point  thus  condensed 
out,  Felippa  [41  ]  obtained  a  set  of  interpolation  functions  for  the  LCCT-12  element. 
These  define  a  piecewise  cubic  polynomial  interpolation  such  that  the  in-plane 
displacements  and  their  first  derivatives  are  continuous  both  in  the  interior  of  the 
element  and  along  the  entire  boundary  ol  the  complete  triangular  element.  A  more 
iletailed  derivation  procedure  and  the  complete  listing  of  the  cubic  interpolation 
functions  are  given  in  Appendix. 

4.2.2  LCCT-11  and  LCCT-9  Elements 

Assembly  of  three  subtriangles  results  in  the  LCCT-12  element  (Figure  3  a). 
However,  the  mid-side  nodal  points  in  this  element  are  not  desirable  for  programming. 
They  complicate  the  mesh  generation  procedure,  increase  the  band-width  of  the 
assembled  equation  systems,  and  require  special  identification  in  calculation  of  the 
stiffness  matrix.  To  overcome  these  difficulties,  it  may  be  desirable  to  develop  a 
special  element  without  external  midpoints.  This  can  be  accomplished  by  assuming  the 
normal  slope  to  vary  linearlv  along  one  or  more  sides  [42]. 

With  the  elimination  of  one  mid-side  node,  the  five-node  element  is  designated  as 
I -(XT- 11  (Figure  3b).  Further  imposing  linear  slope  variation  constraints  on  three  sides 
gives  a  triangle  with  three  nodal  points  and  results  in  LCCT-9  element  as  illustrated 
in  Figure  3(c).  The  LCCT-9  element  is  identical  to  the  Tocher  and  Hartz  [40]  element. 


4.2.3  Quadrilateral  Elements 


Elements  of  quadrilateral  shape  can  be  set  up  as  assemblages  of  triangular 
elements.  Figure  (4)  shows  quadrilateral  elements  built  up  from  four  LCCT-1 2. 
LCCT-1 1  and  LCCT-9  elements.  The  quadrilateral  element  in  Figure  4(a)  has  a  total 
of  23  degrees  of  freedom  for  each  variable  and  was  designated  by  Felippa  as  Q-23. 
Using  four  LCCT-1 1  or  LCCT-9  triangles,  the  Q-19  and  Q-15  elements  as  shown  in 
Figures  4(b).  4(c)  respectix el v  are  realised. 

File  quadnlater.il  element  has  interior  nodal  points  not  connected  to  the  other 
quadrilateral  element  in  a  finite  element  mesh.  These  points  can  be  eliminated  through 
a  local  condensation  process.  Thus,  the  final  quadrilateral  clement  has  24  degrees  oi 
freedom,  corresponding  to  the  two  in-plane  displacement  components  and  their  first 
derivatives  with  respect  to  spatial  coordinates  x  and  y  at  the  four  corners  of  the 
elements  and  an  -*dditional  eight  degrees  of  freedom  corresponding  to  the  normal 
derivatives  of  each  of  the  displacement  components  at  the  mid-side  nodes.  Assuming 
that  the  normal  derivatives  vary  linearly  along  the  edges  of  the  quadrilateral,  the 
mid-side  nodes  can  be  dropped.  This  reduces  the  Q-23  to  Felippa's  Q-19  element  with 
12  degrees  of  freedom  for  each  of  the  displacement  components.  It  is  a  fully 
compatible  quadrilateral  element,  having  a  continuous  cubic  variation  of  displacement 
and  quadratic  variation  of  strain  both  in  the  interior  of  the  element  and  along  the 
entire  boundary  of  the  element,  as  well  as  a  linear  variation  of  normal  slope  along  all 
external  edges.  We  note  however  that  LCCT-9  and  LCCT-1 1  do  not  use  a  complete 
cubic  polynomial.  For  this  reason,  Felippa's  LCCT-1 2  element  based  on  complete  cubic 
interpolation  was  considered  an  improvement  upon  Tocher's  [40]  LCCT-9. 


38 


4.3  Application  to  the  Free-Edge  Stress  Problem 


l;igure  (5)  shows  a  symmetric  laminated  composite  coupon  under  a  state  of 
uniform  axial  strain.  !n  this  case,  away  from  the  ends,  the  transverse  x=constant 
plane  displacement  fields  can  be  assumed  to  be  independent  of  x.  These  assumptions 

imply  the  following  form  for  the  three  components  of  displacement  [  1 1  ]. 
u(x,v,z)  =  ex  +  U(v,z) 

v(x.v.z)  =  \'(v.z)  (75) 

w(x,v.z)  =  \V(  v.z) 

where  e  is  the  uniform  in  plane  strain  in  the  x -direction  and  u,  v.  w  are  components 
of  displacement  along  x,  y  and  /  axes  respectively. 
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4.3.1  Finite  Element  Formulation 


The  constitutive  relationship  for  linear  elastic  anisotropic  material  obeys  the 
generalized  lkxike's  law 

cr  =C  €  i,j=l,2 . 6  (76) 

where  e,  is  namely  the  uniform  extensional  strain  eD.  Based  upon  the  minimum 
potential  energy  principle  (50)  and  substituting  various  interpolation  functions  for 
displacement,  body  force  and  traction  fields  appearing  in  the  governing  functional  (60). 
a  nodal  force-displacement  relation  within  each  element  is  expressed  as 

k  u  =  R  (77) 

'll  I 

where  K,  represents  the  resultant  external  nodal  force.  kn  is  the  element  stiffness 
matrix  which  can  be  written  as 


B  C  B  dV 

im  inn  nj 


m,n  =  l,2, . 6 


The  range  of  i,  j  depends  upon  the  ’degree  of  freedom’  of  the  element,  B  is  the 
displacement  transformation  matrix  and  V  is  the  domain  of  the  element. 

Because  of  the  longitudinal  extensional  strain  is  specified  as  constant,  the 
corresponding  term  in  the  stiffness  matrix  can  be  separated  from  the  rest  and  (78) 
rewritten  as 


k  u  =  k  -R" 

l|  I  I  I 

where  the  range  of  summation  on  m,  n  is  now-  2,  3, 
residual  force  due  to  uniform  in-plane  strain  e„,  i.e. 

R°=  f  B  C  ,edV 

i  I  i  in  ni  1  o 

J  V 


(79) 

6  and  R°  is  the  element 


(80) 
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Alter  forming  t Ire  system  stiffness  matrix  and  nodal  force  vectors,  the  displacement 
components  can  be  obtained  by  solving  the  resulting  set  of  linear  equations  in  the 
standard  manner. 

4.3.2  Higher  Order  Elements 

For  the  free-edge  stress  problem,  due  to  the  fact  that  dependence  of  the 
longitudinal  displacement  on  the  longitudinal  coordinate  x  is  made  explicit,  the  three 
displacement  components  are  completely  defined  by  three  functions  of  two  independent 
transverse  coordinates  y  and  z  as  shown  in  (75>-  Thus,  the  compatible  cubn 
interpolation  functions  used  for  plane  elasticity  problems  can  be  extended  to  the  pseudo 
two-dimensional  model  of  a  laminate  coupon,  and  a  continuous  strain  field  along  both 
in-plane  and  transverse  directions  ensured. 

Figure  (6)  shows  the  nodal  displacement  degrees  of  freedom  considered  in  the 
stiffness  matrix  for  the  complete  triangular  element.  These  included  the  values  of  the 
in-plane  displacement  components  u,  v,  the  transverse  displacement  w,  along  with  the 
first  derivatives  uy,  u2,  v  ,  v2,  w .,  wz  about  the  y  and  z  axes  at  each  corners  i=l,2,3 
as  well  as  the  normal  slopes  at  the  three  mid-side  nodes.  viz. 

U-V  V  un,/  vn4<  vns>  vn(-  Wn,-  Wn5>  This  is  the  I  .(XT- 12  element  but  with  total 

of  36  degrees  of  freedom.  Further  assuming  the  normal  slope  to  vary  linearly  along 
one  or  all  three  sides,  the  l.CCT-11  and  LCCT  9  elements,  with  total  of  33  and  27 
degrees  of  freedom  respectively,  are  obtained  as  specializations. 

As  described  for  plane  elasticity,  quadrilateral  elements,  designated  as  Q-23,  (T19 
and  Q-15,  were  set  up  as  assemblage  of  four  I.CCT-12,  LCCT-11  and  I.CC7T-9 

respectively.  After  eliminating  the  interior  nodal  points  through  a  local  condensation 
process,  the  final  quadrilateral  element,  (J  23,  had  36  degrees  of  freedom,  corresponding 
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to  the  three  displacement  components  and  their  first  derivatives  with  respect  to  the 
spatial  coordinates  v  and  z  at  the  four  corners  of  the  element  and  an  additional  12 
degrees  of  freedom  corresponding  to  the  normal  derivative  of  each  of  the  displacement 
components  at  the  mid-side  nodes.  It  is  a  fully  compatible  quadrilateral  element, 
having  a  continuous  cubic  variation  of  displacement  and  quadratic  variation  of  strain 
not  only  within  the  elements  as  well  as  along  element  boundaries,  but  also  across 
laminate  interfaces. 
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SECTION  V 


CONTINUOUS  TRACTION  FINITE  ELEMENT 
PROCEDURE  FOR  COMPOSITE  LAMINATES 

5.1  Introduction 

In  the  continuous  strain  Q-23  element  developed  for  the  free-edge  stress  problem, 
both  displacement  and  strain  are  continuous  along  interelement  as  well  us  interlaminar 
boundaries.  However,  the  tractions  calculated  across  the  interfaces  between  differently 
oriented  layers  are  discontinuous  due  to  different  orientation  of  adjacent  plies.  Also, 
traction-free  boundary  conditions  associated  with  the  finite-width  laminate  coupon 
cannot  be  satisfied.  In  order  to  remedy  these  two  defects  and  make  the  numerical 
model  more  representative  of  the  real  situation,  it  was  necessary  to  ensure  interelement 
as  well  as  interlaminar  continuities  of  tractions  and  the  traction-free  boundary 
condition  along  the  free-edge.  To  accomplish  this,  nodal  point  degrees  of  freedom  must 
include  some  components  of  stress  and  exclude  normal  gradients  of  displacement  which 
will  be  different  across  interelement  boundaries.  This  was  implemented  by 
transforming  the  displacements  and  their  normal  gradients  at  each  of  the  nodal  points 
of  the  Q-23  element  to  a  mixed  set  of  degrees  of  freedom  w;hich  would  be  continuous 
across  interelement  boundaries.  These  included  both  displacement  and  interlaminar 
traction  components.  Appropriate  displacement-stress  relationships  derived  from  the 
constitutive  laws  were  used.  For  this  element,  traction-free  boundary  condition  could 
be  satisfied  in  a  point-wise  sense. 
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The  continuous  traction  Q  23  element  still  has  cubic  variation  of  displacement  over 
the  element  and  retains  continuity  of  displacements  across  interelement  as  well  as 

interlaminar  boundaries.  The  strains  as  well  as  stresses  vary  quadratically  within  each 

subtriangle  of  the  constituent  LCCT-12  elements  of  the  quadrilateral.  However,  certain 
components  of  strain  are  not  continuous  across  interelement  boundaries  but  interelement 
tractions  are.  This  correctly  allows  for  possible  differences  in  orientation  of  adjacent 

laminae.  If  adjacent  layers  have  the  same  stress  strain  relationships  due  to  identical 

orientation,  stress  continuity  will  imply  strain  continuity  as  well. 

5.2  Derivation  of  Displacement-Stress  Transformation 

In  the  Q-23  element  analysis,  the  number  of  nodal  degrees  of  freedom  is  dilierent 
for  the  corner  nodal  points  and  the  midside  nodes.  For  this  reason,  derivation  of 
transformation  matrices  for  displacements  and  their  gradients  at  the  corners  of  the 
LCCT-12  element,  and  for  normal  ;radients  of  displacement  at  the  midside  nodes  on  the 
element  boundaries,  is  discussed  separately  in  the  following  sections. 

5.2.1  Corner  Nodes 

The  strain-stress  relationship  for  an  orthotropic  material  expressed  in  the  x-y  v 
coordinate  system  is 

ei  =  S|jcri  i.j=l,2,....6  (81) 

or  in  matrix  form 
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es  S,  S,,  S,,  o  O  Slt>  a 

S  s-  S,,  O  o  5,h  °\ 

6z  S,3  *»  S33  °  °  *3t 

=  _  _  182) 
y  0  0  0  S^,  S,.  0  t 

•  yz  44  4>  yz 

y  0  0  0  5tt  0  r 

'  xz  45  55  \z 

y  S.,  S,,  0  0  S,.  t 

•  xy  16  26  26  66  xy 

where  StJ  are  components  of  the  compliance  matrix  for  monoclinic  materials  which 
have  symmem  wnh  respect  to  \  v  plane  0'igure  2)  and  are  defined  as  [4s] 

S, .  =S  m  ’+( 2S,  +S  )nrn'+S,,n4 

I  I  11  12  m*  22 

5  I  2=(  S I  1  +SJ  ~2SI  ~Si„, >m  '  n  ‘  +  S  ,  s 


Sl.t  =  Sl.tm  +SJ3n 


Slfc=[2SI  jin'— 2S2,n‘+(2S(  ,+S(itiKn‘— m')]mn 


S2,=Sj  Jna+(2Si  ,+S(  (  )m^n‘+S,,m4 


^23  ^t3n  +^23m 


§2(i=[2snni-2s22m2+(2S12+S((6xm-~n2)]mn 


S  =S 

^33  °33 

§3fc  =  2(S13_S23)mn 

S44  =  S44m2  +  S5Sn2 
S  =— S.  mn+S  .mn 

4  6  44  6  6 

Vs^y+S^m2 


SM=4(Sl1+S.v“2S.?-SsJm'n2  +  S<, 


where  m=cos0”,  n=sin0’,  and 
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_J_  s 

__ 

"e. 

s  =— 

"ll 

E,,’ 

1-2,  ’ 

13 

E  3  3 

V12 

S  = 

1 

S  =— 

*3  2 

El,’ 

^22 

E22 

^23 

E33 

_^3_ 

s  = 

1 

1  ^ 

In* 

u> 

,  s  = 

1 

E  ’ 
En 

J32 

E22 

’  °33 

E33 

C  —  1  c  —  j  c  —  j 

44  G  ’  55  G  66  G 

23  13  12 

O'  is  ungle  of  ply  from  global  axis  x  to  material  axis  1,  and  Hkk,  Gu,  uit  (i.j.k  =  1.2,3) 
are  moduli  of  elasticity,  shear  moduli,  Poisson's  ratios,  respectively,  in  material 
coordinates. 

Replacing  the  strain  components  by  their  corresponding  displacement  gradients,  lor 
small  strain  theory,  and  rearranging  the  constitutive  relation,  (82)  becomes 

\  5„  5,2  5lb  S,,  0  0  cr, 

\  5,2  5.2  526  S23  0  0  <ry 

UV  ^16  ^26  ^66  S36  0  0  Vv 


Wy  +  VZ 


S13  S23  S3fc  S33  0  0  az 

0  0  0  0  s44  S45  ryt 

0  0  0  0  S45  S5s  rM 


or  svmbolicallv 


€,  =  D.«  D,2  CT1 

e2  n2I  D,,  o\ 


where 


{e  }=  v  {e  }=  w  +v 

y,  i«2/  y  z 


\(T.}=  <rv  ,  tv 


•19 


|€>{€2HD2I]  [Dj  ' 

'{6,} 

(90) 

[s>[d2,Hd2j]  [d„] 

‘[d12] 

(91) 

The  inverse  of  [D,,],  namely,  the  compliance  matrix  in  plane  stress  case,  can  be  'written 
explicitly  as 

o,,  oI2olb 

[Dj  '=  Q|2  Q22  U,b  (92) 

o,„  <3*.  o„, 

where 


5o 


(93) 


0, ,  =0,  ,m4+2(Q12+2QW))m  V+Q22n4 
Qi^=CQ,  ,+0,3— 4Qt>(  >>rn2n2+Ol2(,r‘4+n‘4) 
Q1^=-mn3Q2,+m3nQ1  -mn(m2- n‘XQ12+2Q6fc) 

022=Q,  in4  +  2(0i2  +  2Q66)m2n2+Q22m4 

026=— m3nQ22+mn301  ,+mn(m'— n2XQ12+2Qfc6) 

066=(On+Q22-2012)m2n2+Q6b(m2-n2)2 

with 


'•’ll  1 ' 7  "> 

Q,  J  =  T - LL— .  Q„=T — - - 

^~V  2\V  \2  "  1“,/.2*,2, 


12  21 


Substituting  (92)  into  (91)  and  (90),  (89)  could  be  expressed  as 


w  — B,u  —  B,v  — B.u 

IS,,— X  0  0 

CT 

z  1  x  2  y  3  y 

33 

2 

W  +V 

0  s.„  s,< 

r 

y  z 

44  45 

yz 

U 

0  S,.  S4. 

r 

z 

45  5s 

xz 

where 


^l3^U  +  ^23^12  +  ^3sQl6 
B2_^l3^12  +  ^23^22  +  ^3(.^2(, 


and 


( 04 ) 


(95) 


(96) 

(97) 

(98) 

(99) 


^~B|^13'*’B2^23  +  B3^36 

The  gradients  of  displacement  appearing  in  the  expression  for  interlaminar  strains  were 
then  written  in  terms  of  interlaminar  stresses  using  (95),  i.e. 

(100) 
(101) 


u  =S  AJ  +S,.r 

z  45  yz  55  \z 


V  =(w  +V  )~ W  =S,,T  +S..T  — W 

L  V  /  V  44  V/.  \/  V 
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(102) 


w  =(\\  — 15  u  — B  v  — 15  u  )+15,u  +  15  +I5,ii 
=(S,  —  X)cr  +B,u  +B,v  +15,u 

Combining  (100M102)  with  the  rest  of  the  displacement  nodal  degrees  of  freedom,  u, 
v,  w,  uy,  vy,  wy,  and  noting  that  uN=e0,  the  applied  strain  loading,  the  generalized 
displacement  components  at  the  corner  nodes  of  the  l.CCT-12  element  w'ere  related  to  a 
mixed  set  of  degrees  of  freedom  as  follows: 

O  0  0  0  0  0  O  0 

1  0  0  0  0  O  0  0 

0  1  O  0  0  0  0  0 

0  0  1  0  0  0  0  0 

0  0  0  1  0  0  0  0 

0  0  0  0  1  0  0  0 

0  0  0  0  0  SA.  0 

0  0  0  0  -1  s„,  5,.  o 

4.*>  44 

0  0  B3  B2  0  0  0  S33-X 

or  symbolically 

{r}=[G]{r'}+{R}  (104) 

B,,  B2,  B3,  and  X  occurring  in  (103)  have  been  defined  by  (96)  through  (99).  Thus,  at 
each  of  the  three  corners  in  the  l.CCT  12  element,  we  have  three  displacement 

components  u.  v,  w  along  with  three  inplane  strain  components,  u  ,  v„,  wy,  and  three 

interlaminar  stress  components  rx/,  rV7,  cr7. 
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5.2.2  Mid-side  Nodes 


three  traction  components  which  indeed  can  be  regarded  as  the  out-o!  -plane  shearing 
stress  crnx,  the  normal  stress  crnn,  and  the  in-plane  shearing  stress  crn.  The  relation 
between  xt  and  x',  coordinates  as  shown  in  Figure  (7)  yields 


1»,  =  1,  i12=o,  i13=o 

1 2 1  0»  ^22~ m’  ^23~n 

121— 0,  132  n,  I33— m 


(110) 


n  =0.  n,=m.  11  =11 


where  m=cosO,  n=sin0  and  0  is  the  angle  between  these  two  coordinates.  Substitutin',' 
the  above  quantities  into  (109),  we  have 


t' ,=cr  =mr  +nr 

1  n\  w  \z 


t'  =cr  =nTo-  +2mnr  +n'cr 

2  nn  y  yz  z 


t'  =cr  =(m"~ n  )r  +mn(<T  —cr  ) 

3  ns  yz  z  y 


Here,  for  interlaminar  stresses  acting  on  the  element  boundary  are  nothing  but  the 
traction  field,  crnN,  c rnn,  crns  expressed  in  (112M114)  indeed  represent  the  three  traction 
components  as  illustrated  in  Figure  (7). 

The  stress-strain  relationship  for  the  orthotropic  lamina  expressed  in  the  \-y-x 
laminate  coordinate  system  is 
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\x  here 
and  is 


where, 
axis  1, 


with 


('  is  the  stillness  matrix  lor  monoclinic  symmetry  with  respect  to  x-y  plane 
defined  as  [43] 

CM=C ',,m4+2((:,,+2t ^.Jm^n’+C ^n4 
C 1 2=(C , ,  +C22-4Cfc6)m  2n 2  +C , m  4 + n4 ) 

Ci3-C13m  +C,3n 

cJ?=(:1,n4+2((:12+2rw>)m-,n-+r,2m4 

‘■(rrn 

('  =—  (  mn'+C  m'n—  •<  +2('  hinhn  -n  >  (116! 

it,  . .  ii  i  •  i  • 

(’  =  m  ’n+( ' f  (mn  +K  ^  +2(  ( (  >mnlm'  —  n  > 

r:44=r44ni'+rssnJ 

r:4.s=(CSS-(:44)mn 

f!sS=CS5m2+C44n2 

C66=(C  11  +C  2  2_2C  1  2)n)2 11 2  +C66(  ™  ^  ^ 

with  m=cos0\  n=sin0',  again  0‘  is  angle  between  the  global  axis  x  to  material 
and 


(l — v  v  )F  iv  -¥v  v  )li 

p  _  2.1  32' II  p  _  21  41  2V 'I 

1  I  A  1 2  A 


c.  = 


('.M  = 


A 

l'~*M*, 

,*.V 

A 

(‘]~*i2*2 

,*13 

(:n  = 


( l>  v  v  )l 

1 1  2 1  II  II 


(117) 


r2  1=- 


A  '  -M  A 

^4=G23-  (V=(i,3’  C6(>=Gi2 


A  =  '  ~V  |  2*2 .  ~V 23*32-*}  I  *  1 1  2* 2  I  *3 2*  1 3 


(118) 


(112X114)  along  with  (115),  gives 
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nit  2(>  3<>  22  23  44  4“>  44 
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(119)  is  the  general  relationship  between  traction  components  at  the  element  mill  side 
nodes  and  the  corresponding  displacement  gradients.  It  is  different  for  each  element 
midside  node. 

In  order  to  relate  the  displacement  gradients  at  the  mid-point  to  the  mixed  nodal 
degrees  of  freedom  at  the  two  ends  defining  the  element  surface,  let  j  and  k  denote 
the  first  and  second  cyclic  permutations  of  i= 1,2,3  (i.e.  j=2,3,l  and  k=3,l,2),  the 
projected  dimensions  and  the  corresponding  boundary  length  are  defined  as  (see 
Appendix  A) 

a=yk—  V|(  b=/-z(.  1  =~>J  a'+b'  (120) 

Also,  if  the  outward  normal  is  defined  as  positive  (figure  8).  the  relation  between 
local  and  global  Cartesian  Coordinates  is  [42] 

(si,  a  — b  y— v  ] 
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J 

Considering  ( 1 2 1 )  and  using  chain  rule  of  differentation,  we  have 

=(«H L)  +  =  -2iu 

>.+r  v  Qy  0s,  Qy  Qn  Qy  1  V.r  1  "i+r 
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„  =(4l)  =5“±l  +  -i„ 

z.+3  ve/yi+,  as,  a/  an,  az  1  ■*—  1  " 


1 123) 


1,  Si  +  3  1. 

where  us  3  (i=l,2,3)  denotes  the  tangential  derivatives  of  u  at  one  of  the  mid  side 

nodes  (Figure  9),  which  can  be  further  interpolated  from  the  displacements  and  their 
gradients  at  the  two  ends  of  the  corresponding  element  boundary,  i.e. 


t  a  b  -i  a  b 

u  — u  — — -U  + — -U  +  -t—  u  —  — -u  + — -u 

s,+r  21  1  41  V,  41  /,  21  k  41  >'k  41  zk 


Recalling  (100),  we  havt 


u  —  S -  t  +S,.r 

7.  *'*'  \/  4  -•  N  7 


U.  =S  J .  +S  ,  T 


(124) 


( 125) 

1 1  2o ) 


Substituting  (124),  (125)  and  (126)  into  (122)  and  (1230,  we  can  express  u,,  u,  at  each 
of  the  element  mid-side  nodes  in  terms  of  their  corresponding  normal  derivatives  un,  as 
well  as  the  mixed  nodal  degrees  of  freedom  at  the  two  ends  defining  the  segment. 
That  is 

3a  a'  a  b  a  b  _ 

u  Lu - L_u  +-2-2-S, -T  T  — ) — -  S  ,T 

•vi+3  -m2  i  412  V,  A ,  2  Si.  xz,  A ,  2  4S  v/, 
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yu.+— u  S„r  —  — • —  S  ,r  u 

21“  k  412  >k  41 2  S5  XZk  41 2  4  ”k  1,  n'+a 


(128) 


Following  the  same  procedure  and  considering  (103),  the  transformation  equations  for 
vy,  v,,  wv,  w.  at  each  of  the  element  mid  side  nodes  are 


so 
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o  :  Numbering  for  Mid-Point  Transformation 
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Figure  9:  Numbering  of  Nodal  Points  in  ECCT-12  Element 
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Substituting  (1271(132)  into  (119),  the  relation  between  the  surface  traction  components 
at  each  of  the  element  mid-side  nodes  and  the  corresponding  displacement  normal 
gradients  is 
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crn»  =[T,],vn  +[T2]){r’}i+3  +  {T3}i  (133) 
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md  [Tj,  is  the  product  of  the  transformation  matrix  shown  in  (119) 
following  matrix  [TA ] 
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Rearranging  (133),  the  displacement  stress  translonnatons  lor  the  normal  displacement 
gradients  at  each  of  the  mid  side  mules  in  the  I  ( '( T  12  element  an 


(139) 
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or  symbolically, 

<r„>1+3  =  [^K>,+3+[U<r'}i+3-HS}1 

where 


(140) 


[ii]  -[i  ,l 

1 

(141 ) 

[i.]  =  -[•!•,: 

1,  Mi.] 

(142) 

IS}  =-[T, 

]  M  r,}, 

(143) 

Here,  [ll],  denotes  the  transformation  matrix  which  directly  related  the  normal  slope 
quantities  at  the  mid  point,  i+3,  to  the  corresponding  surface  traction  components.  [L], 
can  be  regarded  as  the  coupling  matrix  between  the  normal  gradients  and  the  mixed 
degrees  of  freedom  associated  with  the  nodal  points,  j  and  k,  at  the  two  ends  of  the 
boundary  {S},  is  the  local  effect  resulting  from  the  applied  uniform  loading  e0. 


5.3  Finite  Element  Formulation 

The  discretized  form  of  (73)  can  be  rewritten  as 
M 

n  =  £(I(q,")T[K,"](qmMqm)TU:"‘})  (144) 

ni-  1 

where  [Km]  is  the  element  stiffness  matrix,  {Fm}  denotes  the  nodal  force  including  the 
local  effect  due  to  the  uniform  in-plane  strain  loading  as  shown  in  (80)  for  the 
free-edge  stress  problem.  {qm}  is  the  set  of  generalized  nodal  displacement  components 
within  the  element  m,  and  M  denotes  the  total  number  of  elements. 
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In  order  to  ha\e  both  displacement  and  traction  continuity  along  interelement  as 
well  as  interlaminar  boundaries,  assuming  (144)  is  based  on  continuous  strain  cubic 
displacement  interpolation  as  in  l.CCT-12  described  in  4.3.2,  the  displacement-stress 
transformation  matrices  derived  in  the  previous  section  were  imposed  on  the  generalized 
degrees  of  freedom.  Combining  (104)  and  (140),  the  generalized  displacement 
components  in  the  l.CCT-12  element  were  transformed  to  the  mixed  type  degrees  of 


freedom  as  Id  lows: 


I'M 


!4<l» 


w  here 


[Tm] 


[ti]  n  o  o  o  o 

o  [(,]  o  o  o  o 

0  0  [(»]  0  0  0 

iu3  uj3  o  im3  o  o 

o  [L],  lU,  o  [H],  o 

[1J2  o  UJ,  0  0  [ll]2 


(145) 


(146) 


(147) 


<q,',i  =  {u1,v1,wI,uV),v>i,wVi.u/|,v/i.wVi,u2.v;).w,u>2.vy2,wy2,uz2,vz2,w^ 


U  ,V  ,W  ,11  A  ,W  ,U  A  ,W  .11  A  AV  ,U  A  AV  ,U  ,V  ,W  I 
'  '  '  d  M  ",  ,:4  "4  "4  "s  "s  "s 


(1  IS) 


l<rv  =  )u,A  ,.W  ,4»S|.vV(.wNi.Ts/i.7Wi.a|.u;,A  ,.w2.ux/vv/wv/T^,Ty^,a  /U(Arw  r 

U  A  ,W  ,7  ,7  ,<T  ,CT  ,a  ,C7  ,CT  ,CT  ,CT  ,CT  ,<7  ,C7  }'  (149) 

V  j  V3  v3  S',3  VZ3  /-J  ||'G  ns,,  "'4  '>"4  ns4  n>“,  n,15  ,lsS 

liach  of  the  [!.],  matrices  in  (146)  was  divided  into  two  parts  to  match  the  mixed 
degrees  of  freedom  in  (149).  In  the  finite  element  computation,  this  transformation 
was  implemented  during  the  formation  of  each  of  the  subtriangular  element  stiffness 


matrix  and  load  \ectoi  corresponding  to  the  I.CCT12  element. 


Thus,  due  to 


appearance  of  displacement  and  interlaminar  traction  components  at  the  corners  of  the 
triangle  as  well  as  traction  components  at  the  mid-side  nodes,  a  cubic  variation  of 
displacement  and  a  quadratic  variation  of  traction  were  ensured  along  element 
boundary.  More  importantly,  both  fields  are  continuous  across  the  common  boundary 
between  two  contiguous  triangular  elements.  Substituting  (145)  into  (144),  we  have 


n  =  £<> 


,ni\Tr.rniiT 


)  [k'"]  ([T",]iq-i+{i>",i)— (urn-i 


r){rr"}) 


=  £( \k)'  I'l  I"  ]‘[K"]  criun  +  llq-'  }‘[  I  I'lh"  ]ll'  "i  +  -^-)l)1"}l|K"‘]  IT'  J{q" 


+iir"')V"](p'"Hq" 


Trr M]y  >",}In 


or 

M 

ft  =  £<  |{q'm}T[Knl]{q'm}-{q'rn}T{Rm}+Cm)  ( 1 51 ) 

m=  1 

w'here 

[K">[T,nf  [K,n]  [Tm]  (152) 

{ H ,n } =[Tm]T{  F"‘}— (T"']Tf  K '"]{ Pm }  (153) 

C1" = ■—  { P"  }T[  K  "’]( P'1}— { P'"  |T{  l'l:i  | =constan  t  (154) 

Using  global  coordinates,  (151)  could  Ire  written  as 

ft  =  j  {q'}T[k]{q'Hq'}{R}+Cni  ( 1 55) 

Taking  the  variation  of  (155)  with  respect  to  {q*}  yields  the  system  generalized  nodal 
force-displacement  equilibrium  equations 

[FTKq'MR}  (156) 
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I'he  displacement  as  well  as  interlaminar  stress  components  were  then  obtained  by 
solving  the  resulting  set  of  linear  equations  in  the  standard  manner. 

5.4  Calculation  of  Stresses 

Solutions  of  the  finite  element  system  gives  the  three  displacement  components, 
their  tangential  gradients  along  the  element  edges,  rotation  about  the  longitudinal  axis 
and  interlaminar  stress  field  at  the  corner  nodes  of  t he  Q  23  element,  along  with  the 
boundary  traction  components  at  the  center  point  oi  each  of  the  Q  23  element  surfaces. 
To  retrieve  the  rest  of  the  displacement  components  at  each  corner  node,  the 
transformation  matrix  used  in  ( 104)  was  reapplied  to  the  calculated  nodal  point 
solution.  Having  found  the  complete  displacement  gradients  and  hence  strains,  the 
inplane  stress  components  at  the  corner  nodes  of  the  Q-23  element  could  then  be 
computed  using  (115). 

The  interlaminar  or  interelement  stress  components  at  the  mid-point  of  each 
element  boundary  are  merely  the  traction  components  directly  produced  by  the  solution 
of  the  continuous  traction  Q-23  element  provided  a  rectangular  mesh  is  used  in  the 
analysis.  For  determining  the  rest  of  the  stress  components,  the  normal  displacement 
gradients  at  mid-side  nodes  were  recovered  from  the  nodal  point  solution  using  (140). 
Furthermore,  the  displacement  and  their  gradients  along  the  edge  at  the  mid-point  w’ere 
interpolated  from  the  previously  computed  displacement  components  at  the  tw'o  ends  of 
the  segment.  Having  transformed  these  displacement  gradients  from  local  to  global 
Cartesian  coordinates,  calculation  of  the  remaining  strain  and  stress  components  at  the 


mid-side  of  each  of  the  element  boundaries  is  direct. 


5.5  Boundary  Conditions  of  a  Quadrant  of  the  Delamination 
Specimen 

For  the  free-edge  delamination  specimen,  because  of  symmetries  in  the  laminate, 

only  one  quadrant  of  an  x=constant  plane  was  considered  (Figure  2).  Along  the 

boundary,  either  displacements  or  tractions  are  specified  at  each  point. 

5.5.1  Boundary  Conditions  Along  Lines  of  Symmetry 

Symmetry  of  loading  as  well  as  geometry  about  the  mid-plane  implies  that  the 

displacement  functions  satisfy  the  following  condition 

f.  (y,-/.)= l  (y, z )  (157) 

Y(\V/J=V(y.z)  ( 1 5S) 

\V(y,-z)=-W(v,z)  (159) 

Using  chain  rule  of  differentiation, 

— U2(y,-z)=U2(y,z)  (160) 

— Vz(y,-z)=Vz(v,7.)  (161) 

Wy(y,-z)=-Wv(y,z)  (162) 

Setting  z— 0  in  (5.79)-(5.82),  we  have 

\V(y,())=0  (163) 

\Y  (v,0)=L  (v,0)=V  (v,0)=0  (164) 

V  *  7.  '  /.  ■ 

From  (5.84),  'yX7(y,0)=y  (y,0)=0,  and  consequently,  for  the  layered  orthotrnpic  material. 

r  (v,0)=r  (y,0)=0  (165) 

xz  -  1  yz 

Also 

(Wv-VzXv,o)  =  o 
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Invariance  under  a  rotation  of  180  degree  about  the  /  axis  through  the  center  of  the 
specimen  implies 


L'(-y,z)=— U(y,z)  (166) 

V(-y,z)=-V(v,z)  (167) 

W(-y,z)=W(y,z)  (168) 

(166)  and  (167)  lead  immediately  to 

L(0,z)=Y(0./)=0  ( 169) 


lor  all  /..  and  consequent  I  v 

l  (0,z)=\  (0./)=0  (170) 

/  / 

By  chain  rule  of  differentiation,  (166)  yields 

-Wy(-y,z)=\Vy(v,?.)  (171) 

Hence,  for  y=0 

Wy(0,z)=0  (172) 

(170)  and  (172)  imply  -yx2(0,z)=yJ2(0,z)=0.  Hence,  for  the  layered  orthotropic  material, 
Txz(0,z)=ryz(0,z)=0  (173) 

Also 

(W  -Y  )(<)./)  =  o 

(-ombining  (163).  (164).  (165),  (169).  (172)  and  (173),  along  with  U(0.0)=0  in  order 
to  present  rigid-body  displacement  of  the  laminate,  the  continuous  traction  finite 

element  model  for  a  quadrant  laminate  under  consideration  should  satisfy  the  following 
conditions  along  lines  of  symmetry 

L’(0,0)=U(0,z)=V(0,z)=W(y,0)=W  (v,(»=\V  (0,z)=0  (174) 

y  -  y 

T  (v,0)=T  (\’,0)=T  (O.Z )  — T  (0,Z)=0 

VZ  ‘  \f  V/ 
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(175) 


5.5.2  Tract  ion -Free  Boundary  Conditions 

The  traction  !  ree  boundary  conditions  associated  with  tine  quadrant  of  the  laminated 
specimen  are 

r  (\\H)=t  (v.H)=o- (y,H)=0  (176) 

at  the  top  surface,  along  with 

cr  (B,/.)=r  (H,z)=0  (177) 

\  ’  V  V  W 


at  the  lateral  1  lee  edge.  Ilete  2H  and  211  denote  the  total  width  and  thickness 
respite 1 1\ el v  ol  the  laminated  specimen.  1  smg  the  lontmuous  traction  (.)  23  model. 


tract  ion -free  b.uindaiw  conditions  shown  m  <  1 7<>  ■  ian  be  identicaiiv  satisfied  for  nodal 
points  on  the  top  surface.  However.  due  to  the  lack  of  inplane  stress  components  as 
nodal  degrees  of  freedom,  only  the  last  condition  shown  in  (177)  can  be  specified  at 
those  element  corner  nodes  along  the  lateral  free-edge.  To  completely  satisfy  the 
traction-free  condition  along  the  lateral  free-edge,  the  following  device  was  developed. 

In  older  to  enforce  the  remaining  two  inplane  stress-free  conditions  in  (177),  it  is 
necessary  to  express  these  in  terms  of  nodal  point  degrees  of  freedom.  This  results  in 
a  linear  relationship  Iretween  the  degrees  of  freedom  at  nodal  points  on  the  free  edge 
1'he  stress  f-ee  condition  can  be  written  explicitly  as 

cr  =()=('  e  -K\.w  +(',,\  u  (178' 

t  =()=(’  e  +(\  w  +C  u  <170) 

or  symbolically 


M 

c  c  c  c 

12  23  22  26 

"1 

c  c  c  r 

If,  36  26  :*■ 

(180) 
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where  ('.  <ire  components  of  the  silliness  matrix  del  ined  by  (116). 
ecjuation.  the  inplane  strain  components  n,  %.  can  be  expressed 
interlaminar  normal  strain  w„  by 
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Solving  above 
in  terms  of 


(181) 


w  here 


1  =!(('• ,  ( ‘  -('  (',  > 
a 


( 1 8  J  » 


and 


Q  ^'22^'bb  ^2b 


(183) 


Substituting  (102)  into  (181)  for  w7,  we  can  further  relate  uy  and  v  to  the 
interlaminar  normal  stress  component  cr7  through  the  following  linear  relationships 

u  =p,cr  +p7e  (184) 

S'  •  1  1  I  2  o 


\  .  =q  cr  .+<)  ,e 


(185) 


w  here 


I’rjlCA  -vl 

P,2'j<lnl22Bi-|,;l..,B;+l..2B,  +  l2,) 

li'jl'iAj-X'I 


(186) 

(187) 

( 1 88) 

(189) 
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and 


(190) 


0=1~1,,» 

Again,  1),,  11,,  15^  and  \  occurring  in  (186)  through  (190)  have  been  defined  by 
(96)-(99).  Thus,  the  variables  associated  with  the  nodal  points  of  the  free-edge  in  the 
lateral  continuous  traction  Q- 23  element  are  no  longer  independent  but  related  through 
(184)  and  (185).  Incorporating  these  linear  relationships  in  the  displacement-stress 
transformations  shown  in  (10-4)  and  (14(0  lor  elements  on  the  lateral  boundary  implies 
satisfaction  of  the  traction  live  boundary  conditions  (177).  I  he  transformation  becomes 
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where 

r,  =  S,,— X+B^+B^p, 
T ,  ~  B1  +  B,q,+ll)p, 


(192) 

(193) 
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SECTION  VI 


ANALYSIS  OF  FREE-EDGE  DELAMINATION  IN 
LAMINATE  COMPOSITE  SPECIMENS 

6.1  Introduction 

Continuous  traction  finite  element  formulation  developed  in  the  previous  section 
was  applied  to  obtain  an  approximation  to  displacement  as  well  as  stress  fields  in  a 
free-edge  delamination  specimen.  The  analysis  consisted  of  two  parts.  The  first 
consisted  of  solving  four-ply  symmetric  laminates.  The  purpose  was  to  examine  the 
credibility  of  the  continuous  traction  finite  element  model  by  comparing  the  numerical 
solutions  with  those  from  Pagano's  analysis  based  on  a  generalization  of  Reissner’s 
theory.  At  the  same  time,  the  deficiency  in  using  the  continuous  strain  free-edge 
stress  model  described  in  Section  IV  was  examined.  The  second  part  consisted  of 
studying  of  edge  delamination  tendency  in  two  classes  of  multi-ply  laminates  subjected 
to  longitudinal  loading.  Tor  laminate  specimens  with  stacking  sequence  of 
[(0/— 0)m/9On  jl.,  displacement  and  stress  fields  calculated  from  the  continuous  traction 
Q-23  element  were  compared  with  those  obtained  using  an  overlay  procedure  with 
constant  strain  elements  [45]  and  with  experimental  observation  [46].  For  the 
[(0/— 6)3/90],  laminate  specimens,  the  continuous  traction  finite  element  code  was  used 
in  conjunction  with  some  well  known  anisotropic  failure  criteria  [47-50]  to  predict  the 
onset  of  transverse  cracking  and  the  onset  of  edge  delamination  in  various  laminate 
specimens  observed  in  experimental  data  [51].  The  purpose  was  to  evaluate  the 
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suitability  of  these  failure  criteria  lor  application  to  the  free-edge  delamination  problem. 
Because  ol  symmetries  in  the  laminates,  only  one  quadrant  (Figure  2)  was  considered 
in  each  case. 


6.2  Four-Ply  Laminates 

In  this  section,  analysis  of  two  long  symmetric  laminate  strips  made  of 
graphite  epoxy  materials,  with  liber  orientations  of  [-15  45]  and  [090]  under  uniform 
inplane  strain  in  the  longitudinal  direction  is  described.  The  relation  between  laminate 
width  and  thickness  was  2b  l(>h  lollowmg  [ 8 ].  In  the  analysis,  each  plv  was 
idealized  as  a  homogeneous,  elastic  orthotropic  material.  lor  comparison  purpose,  the 
material  properties  assumed  here  following  Pagano's  work  [8] 

H, ,  =20xl()bpsi 
I:2,=E31=2.1  Xl06psi 
G  j  2 =G  j  3 =G  23  =0.85  X 1 06  psi 

,/12=,'«3=»'23=0-21 

The  subscripts  1,  2  and  3  correspond  to  the  longitudinal  transverse  and  thickness 
directions  respectively.  A  144-element  model  as  shown  in  Figure  (10a)  was  used  to 
discreti/.e  a  typical  \=constant  plane.  Numerical  results  based  on  the  continuous 
traction  (,)  23  and  continuous  strain  (,)  2 5  elements  were  compared  with  Pagano's  [8] 
analytical  solution. 

Figures  (11)  through  (23)  illustrate  comparisons  for  both  stress  and  displacement 
fields  at  specific  locations  for  the  angle-ply  and  cross-ply  laminates  using  different 
solution  techniques.  The  value  of  N  in  these  figures  corresponds  to  the  number  of 
sublayers  used  in  Pagano's  theory.  Thus,  N=6  indicates  that  each  physical  layer  of 
thickness  h  was  modeled  by  three  sublayers  each  of  thickness  h  3,  while  N  =  2  denotes 
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that  each  physical  layer  i.s  treated  as  a  unit  |S].  Also,  the  calculated  displacements  and 
stresses  were  normalized  by  the  applied  unilorm  strain  loading  e,.  which  has  been 

taken  as  unit  in  the  present  analysis. 

6.2.1  Angle-Ply  Laminate 

Figures  ( 1 1 )  and  (12)  show  the  distribution  of  o\  and  along  the  width  of  the 
laminate  at  the  center  line  of  the  top  (45  degree)  layer.  The  results  obtained  iisinc 
the  continuous  traction  (J  25  element  agreed  ipnte  well  with  those  of  Pagano's  \  <> 
solutions  across  the  entire  width  ol  the  laminate. 

.A  comparison  of  the  shear  stress  (rs 4  distribution  along  the  interface  ol  the 

45  -45  layers  (Figure  13).  indicated  that  the  continuous  traction  Q  23  solution  had 
sharp  rise  toward  the  free  edge  similar  to  Pagano's  solution  with  N=6.  Satisfactory 

agreement  was  observed  between  these  two  solutions  for  stress  across  the  width  except 
at  the  free-edge  boundary  where  continuous  traction  Q  23  somewhat  underestimated  the 
singular  stress.  Figure  (14)  show's  the  through-thickness  stress  distribution  of  i\ 
calculated  from  both  continuous  strain  and  continuous  traction  Q-23  elements  at  the 
free-edge  of  the  laminate.  Very  ciose  agreement  was  generally  observed  between  these 
two  solutions  throughout  the  thickness.  Also,  at  the  interface  of  the  45  45  layers, 
continuity  of  the  interlaminar  shear  stress  was  ensured  for  both  cases.  This  is  because 
of  the  rotational  symmetry  about  /  axis  in  the  particular  angle-plv  laminate  considered. 
The  singular  behavior  of  which  is  highly  localized  at  the  interface  between  45  45 

layers,  is  noticeable.  The  distribution  of  oy  along  the  interface  between  45°  and  -45" 
plies,  which  was  not  indicated  in  (8],  will  be  discussed  in  the  next  section. 

For  the  axial  displacement  distribution  across  the  width  of  the  top  surface, 

continuous  traction  (,)  23  results  compared  well  with  Pagano's  N  -(>  solution  (figure  15). 
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CONTINUOUS  TRACTION 
PfiGRNQ  (N  =  6) 


Figure  (16)  shows  the  through-thickness  distribution  of  axial  displacement  based  on 
continuous  strain  as  well  as  continuous  traction  Q  23  elements  at  the  free-edge  of  the 
laminate.  Again,  these  two  solutions  matched  well  throughout. 

6.2.2  Cross-Ply  Laminate 

Distribution  of  tr.  along  the  width  on  the  central  plane  of  the  [0/9()l.  laminate, 
shown  in  Figure  (17).  indicates  a  sharp  rise  near  the  free -edge  boundary.  Solutions 
obtained  from  the  continuous  traction  Q  23  element  nearly  coincided  with  I'ai'ano's  \  t> 
solution  twer  the  entire  width  ('!'  the  laminate. 


Figure  (18) 

shows  the  variation  of  cr.  along 

ihe 

interface 

between  the 

O'  anil  9(> 

plies.  Due  to 

the  presence  o 

1'  the  discontinuity 

in  e 

lastic  * 

-  'ties. 

a  sin; 

gular  stress 

behavior  would 

l>e  expected 

at  the  f  ree  edge. 

On 

the  a1 

•ary. 

result 

from  the 

continuous  traction  Q-23  element  had  a  steeper  gradient  than  that  of  Pagano’s  theory. 
Apparently,  one  possible  reason  for  this  discrepancy  is  that,  in  Paganos  analysis,  each 
physical  layer  could  be  modeled  by  at  most  three  sub-layers.  However,  in  the  finite 

element  analysis,  the  thickness  was  divided  into  eight  elements.  If  a  coarser 

discretisation  were  to  be  used,  say  4x18,  c ri  calculated  from  the  continuou  traction 
Q-23  element  would  possibly  agree  quite  well  with  Pagano’s  solution  at  the  free  edge 
interface.  Figure  (19)  shows  the  influence  ol  through  t be  thickness  refinement  of 

mesh  on  a  .  Figure  (20)  shows  through  the  thickness  distribution  of  crj  at  the 

f ree-edge  of  the  laminate.  In  the  vicinity  of  the  interface,  the  continuous  traction 
Q-23  element,  enforcing  continuity  of  cry  at  the  interface  between  differently  oriented 
layers,  gave  a  stress  distribution  quite  different  from  that  given  bv  the  continuous 
strain  Q-23  analysis.  Away  .  ,om  the  interface  the  two  sets  of  results  were  close. 
Also,  the  interlaminar  stress  cr  was  observed  to  base  a  maximum  value  in  the  interior 
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of  the  9o  deg  layer  closer  to  the  interlace  with  the  top  layer.  However,  in  botli  cases, 
the  solutions  displayed  oscillatory  patterns  near  the  interlace.  This  could  Ire  due  to  the 
finite  element  mesh  used  being  not  fine  enough  to  approximate  the  steeply  varying 
stresses  associated  with  abrupt  change  of  material  properties. 

Values  of  t  along  the  interface  between  the  [0/901,  layers,  calculated  from  the 

continuous  traction  Q-23  element  (Figure  21).  showed  satisfactory  agreement  with  those 
calculated  by  l'agano.  1'his  is  because  the  continuous  traction  Q-23  element  exactlv 

satisfies  the  traction  free  boundary  condition  similar  to  l’agano's  theory.  However,  an 
oscillatory  error  was  observed  near  the  I  ree  edge.  Apparently,  further  mesh  relinement 
along  the  y-direction  is  required  near  the  free-edge  in  order  to  approximate  the  singular 
stress  behavior.  Figure  (22)  displays  through-the-thickness  stress  distribution  ol  r, 

calculated  from  both  continuous  strain  and  continuous  traction  Q-23  elements  at  the 

free-edge  of  the  laminate.  Apparently,  satisfaction  of  the  traction-free  boundary 
condition  associated  with  the  continuous  traction  Q-23  element  represents  an 

improvement  over  the  continuous  strain  Q-23  element. 

Comparative  results  for  the  variation  of  transverse  displacement  along  the  top 

surface  of  the  [0/901  laminate  are  shown  in  Figure  (23).  Fxcellent  agreement  was 

observed  between  results  using  the  continuous  traction  Q-23  element  and  Pagano's  N-2 

solution. 
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6.2.3  Effect  of  Traction-Free  Edge  on  the  Solutions 

In  order  to  investigate  the  effect  of  requiring  the  satisfaction  of  a  traction-free 
boundary  condition  on  the  finite  element  solutions,  the  continuous  traction  Q-23 
element  was  employed  with  only  the  requirement  that  Tyi=0  along  the  lateral 
free-edge  of  the  four-ply  laminate  specimens.  In  other  words,  the  displacement 
constraint  conditions  developed  in  (188)  and  (189)  used  to  specify  the  in-plane 
stress-free  boundary  conditions  were  not  imposed  in  this  model,  for  convenience  in  the 
following  comparisons,  this  is  designated  as  continuous  traction  (partial). 

figure  (24)  shows  the  distribution  of  r,  along  the  interface  of  the  45/-45  layers. 
Solutions  calculated  from  the  continuous  traction  (partial)  had  a  steeper  gradient  in  the 
vicinity  of  the  free-edge  than  that  of  previous  continuous  traction  Q-23  element.  A 
similar  observation  is  made  for  the  variation  of  tr7  at  the  interface  of  [0/901,  laminate 
(Figure  25).  Thus,  it  is  concluded  that  the  nonimposition  of  the  conditions 
<r  =0  and  t  = 0  would  overestimate  the  magnitudes  of  the  interlaminar  stresses  on  the 
interface  near  the  free-edge.  However,  the  nonsatisfaction  of  these  two  traction-free 
boundary  conditions  had  no  significant  effect  on  the  displacement  field  in  the  laminate 
specimens.  Figure  (26)  shows  the  solutions,  for  axial  displacement  distribution  across 
the  width  of  the  top  surface  in  the  [45  -45].  laminate,  obtained  from  the  continuous 
traction  (partial)  and  from  the  continuous  traction  Q-23  element  which  satisfies  all 
traction  free  conditions  at  the  free-edge. 
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6.2.4  Effect  of  Mesh  Refinement 


The  analysis  of  the  four-ply  laminate  specimens  was  originally  carried  out  by- 
using  the  continuous  traction  Q  23  element  with  a  uniform  80-element  model  as  shown 
in  Figure  (27).  However,  for  reliability  of  results,  the  finite  element  mesh  must  be 
refined  in  regions  of  steeply  varying  stresses.  This  results  in  the  present  analysis 
using  the  144-element  model  which  indeed  was  obtained  by  dividing  the  two  elements 
closest  to  the  free-edge  in  the  80-element  model  into  10  elements  along  the  y -direction. 
To  study  the  effect  of  mesh  refinement  on  the  continuous  traction  finite  element 
solutions,  comparisons  were  made  for  the  stress  distributions  in  the  lour  ply  laminate 
specimens  between  the  uniform  80-element  and  the  locally  refined  144-element  models. 

Figure  (28)  shows  the  distribution  of  rv;  at  the  interface  of  [45 '-45].  laminate 
based  on  the  144-element  model.  A  steeper  gradient  of  rx;  was  observed  on  the 
boundary  as  compared  with  the  result  using  80-element  model.  A  comparison  of  cr7 
distribution  along  the  interface  of  the  45/-45  layers,  indicates  that  the  144-element 
model  had  a  compressive  finite  maximum  value  at  the  free-edge  rather  than  a  tensile 
quantity  from  the  80-element  model  (Figure  29).  This  indeed  has  demonstrated  the 
inappropriate  sign  of  cr;  shown  in  Figure  (l)  based  on  the  perturbation  technique  as 
well  as  finite  difference  method.  For  variation  of  cr;  along  the  interface  of  [0/90] 
laminate.  Figure  (30)  indicates  that  a  singular  stress  behavior  was  properly  reproduced 
by  the  144-element  model  and  did  not  show  well  in  the  results  based  on  the 
80-element  mesh. 

Use  of  144-element  model  might  still  be  insufficient  to  approximate  the  singular 
stress  behavior.  One  example  is  the  ry7  distribution,  which  had  an  oscillatory  pattern 
near  the  free-edge  along  the  interface  of  [0/90\  laminate,  as  mentioned  before.  To 
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overcome  this,  a  more  refined  mesh  (20S  element  model'  obtained  by  further  dividing 
the  two  elements  closest  to  the  free-edge  in  the  144-element  model  into  10  elements, 
was  used  in  the  analysis.  Along  with  the  results  calculated  from  80-  and  144  element 
models  respectively.  Figure  (31)  indicated  the  improvement  of  the  rvz  distribution  along 
0/90  interface  over  the  boundary  layer  region  as  more  refined  elements  were  used  near 
the  free-edge.  However,  regardless  of  mesh  patterns  used,  there  was  an  oscillatory 

error  in  r.  in  the  two  elements  next  to  the  free-edge. 

Figure  (32)  shows  through  the  thickness  distributions  of  a  at  the  free  edge  ol 

lo  90]  laminate  based  on  the  144-element  model  but  with  finer  mesh  neai  the 

interface  (Figure  10b)  and  its  refinement  (288-element  model)  in  the  thickness  direction. 
It  is  observed  that  the  oscillatory  error  near  the  interface  was  reduced  by  using  the 

refined  144-element  model  and  was  nearly  disappeared  under  more  refinement  over  the 
laminate  thickness.  Meanw'hile,  the  maximum  value  of  a,  within  the  90-degree  layer 
was  moving  closer  to  the  0/90  interface. 
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6.3  Free-Edge  Delamination  in  Multi-Ply  Laminate  Specimens 

Analysis  of  the  tour-ply  laminate  specimens  described  in  the  previous  section 
demonstrated  some  validity  of  using  the  proposed  finite  element  procedures  in  solving 
free-edge  effect  problems.  Both  continuous  strain  and  continuous  traction  Q-23  elements 
had  similar  prediction  on  the  displacements  and  inplane  stress  distributions,  which  also 
compared  well  with  Pagano's  analytical  solutions.  However,  discrepancy  between 
continuous  strain  and  continuous  traction  finite  element  models  was  apparent  lor  the 
interlaminar  stresses  near  laminate  interlaces  between  differently  oriented  layers  or  near 
the  traction-lree  boundary.  Hue  to  the  fact  that  stress  continuity  across  interlaminar 
boundary  as  well  as  traction-lree  boundary  condition  are  exactly  satisfied  in  the 
continuous  traction  Q-23  element,  solutions  obtained  from  this  approach  were  expected 
to  be  more  reliable  than  those  from  the  continuous  strain  Q-23  element.  In  this 
section,  application  of  the  continuous  traction  Q-23  element  to  investigate  the  free  edge 
effect  as  well  as  initiation  of  edge  delamination  in  the  multi-ply  laminates  is  described. 

6.3,1  Analysis  of  [(0/-0)m/9Oy2]s  Laminates 

Four  types  of  laminates  with  predetermined  fiber  orientations  [46]  were  used  in 


present  : 

investigation.  These 

are 

Type 

Stacking  Sequence 

Width 

Ply  thickness 

Plies 

A 

[(49.8/-49.8),/90]v 

1.0  in 

0.(K)506  in 

22 

B 

l(30.8/-30.8)s/90]s 

1.0  in 

0.00508  in 

22 

C 

[(25.5/-25.S)s/90]s 

1.0  in 

0.00505  in 

22 

D 

[(47.9/-47.9)JO/90]s 

1.0  in 

0.00499  in 

42 

The  material  used  in  the  study  was  AS4/3501-6,  graphite-epoxv,  and  the  elastic 
constants  were  [46] 
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K( 

1 .32x  lo'psi 

G  J  ,=0.83x1  <)'’psi 
vt  ,=0.35 

All  the  specimens  have  been  investigated  both  analytically  and  experimentally  at 
the  AFWAI./AFFDI.  [46].  The  Delamination  Moment  Coefficients  (DMC)  were  derived 
and  used  to  evaluate  quantitatively  the  delannnation  tendency  of  the  laminates.  Also, 
generalized  constant  strain  element  was  applied  to  analyze  half  of  the  width  of  the 
laminate  specimens.  In  the  experimental  aspect,  xarious  techniques  including  Transxerse 
Strain  Gages,  (’racked  Silver  Ink  Instrumentation  and  Acoustic  Emission  Instrumentation, 
etc.  were  used  to  determine  the  onset  of  delamination  and  to  validate  the  analytical 
results. 


6.3.1. 1  Numerical  Evaluation 

A  154-element  model  shown  in  Figure  33(a)  was  used  to  discretize  one  quadrant 

of  a  typical  x=constant  plane  in  the  laminate  specimen  Types  A,  B,  C,  and  a 

294-element  model  (Figure  33(b))  w'as  used  for  speicmen  Type  D.  Each  ply  was 

modeled  by  a  single  element  through  its  thickness.  Interlaminar  stress  field  within 
various  laminate  specimens  for  an  applied  longitudinal  average  stress  of  1(X)  psi  were 
computed. 

Comparisons  of  c rl  distribution  along  the  mid-plane  of  various  multi-ply  laminates. 
(Figures  34-37)  indicate  that  the  continuous  traction  Q-23  solutions  had  sharp  rise 
toward  the  free-edge  similar  to  constant  strain  element  solutions  [46].  Satisfactory 
agreement  was  generally  observed  between  these  two  solutions  for  stresses  in  the 

vicinity  of  the  free-edge  except  on  the  boundary  where  the  stress  calculated  using  the 
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Q-23  element  was  distinctly  less  than  that  front  the  constant  strain  solution. 
However,  the  cr,  values  calculated  front  the  constant  strain  element  were  extrapolated 
from  the  interlaminar  stresses  at  z=0  obtained  by  l.agrangian  interpolation  of  the  o\ 

values  at  the  element  centroids.  This  is  unlike  the  continuous  traction  Q-23  solution 

where  <x;  at  the  free-edge  was  directly  calculated  as  nodal  degree  of  freedom  and 

would  be  expected  to  be  more  reliable. 

figures  (3SM41)  show  the  through-tile  thickness  stress  distributions  of  cr  and  r, 
calculated  from  continuous  traction  Q  23  element  at  the  free-edge  of  \arious  laminate 
specimens.  It  is  observed  that  for  the  same  applied  axial  stress,  specimen  D  had  the 
largest  value  of  normal  stress  cr,  at  the  tree -edge,  followed  by  specimens  A,  B  and  ('. 
Figure  (42)  illustrates  this.  The  slope  discontinuity  of  cr,  at  the  interfaces  of  the 

free-edge  shown  in  Figures  (38W41)  was  possibly  attributed  to  the  material  as  well  as 
geometrical  discontinuities  in  that  region.  Figure  (43)  indicates  that  much  smoother  a7 
distributions  through  the  laminate  thickness  were  recovered  within  the  angle-ply 
laminae  at  a  small  distance  from  the  free-edge.  In  fact,  with  further  refinement  along 
the  free-edge.  the  solution  of  cr.  on  the  free-edge  is  even  better.  Figures  (44)-(45) 
show  the  solution  for  op  at  the  free-edge  and  at  y=0.495  for  specimen  A  with  each 

edge  element  being  refined  into  four  elements  along  y  direction.  Also,  figure  (46) 

shows  the  functional  dependence  of  longitudinal  stress  o\  as  well  as  interlaminar 

normal  stress  cr;  on  the  fiber  orientation,  respectively,  for  the  [(0/— 0)5/9Ol.  laminate 
under  the  same  applied  loading.  The  ordinates  of  these  curves  are  the  respective 
values  of  crx  and  o\  at  the  intersection  of  the  mid-plane  and  traction-free  edge  of  the 
laminate.  Results  obtained  from  the  continuous  traction  Q-23  element  indicated  that 

troth  cr,  and  <r,  attained  their  maximum  values  approximately  in  the  fiber  orientation 


111 


6  =  30°. 


The  shear  stress  rvi  distributions  were  similar  for  these  specimens,  and  their 
magnitudes  are  relatively  smaileT  than  the  maximum  normal  stress  cr2.  However,  the 
existence  of  rXJ  evidently  reflected  a  defect  of  the  numerical  model  adopted  in  [46]  in 
which  the  delamination  specimen  was  treated  as  an  axially  symmetric  problem  in 
which  rVi  was  inherently  assumed  to  be  zero  throughout  the  laminate  thickness.  It 
was  noted  that  o\  distribution  had  a  slope  discontinuity  at  the  mid  plane  surface 
within  the  90-degree  layer.  This  does  not  appear  to  be  reasonable  lor  the  present 
symmetric  laminate  specimens.  Presuming  that  this  was  associated  with  the  use  of  a 
single  element  through  the  thickness  of  90-degree  layer  being  insufficient  to 
accommodate  the  mismatch  at  the  interface  between  angle-ply  and  cross-ply  laminae,  a 
study  was  carried  out  refining  the  mesh  near  the  midplane.  Figure  (47)  shows  the 
dramatical  improvement  of  cr1  distribution  near  the  mid-plane  surface  of  specimen  A 
as  increasing  number  of  elements  was  used  in  the  discretization  of  90-degree  layer.  At 
the  same  time,  the  maximum  value  of  crt  in  the  interior  of  the  transverse  layer  was 
observed  to  move  closer  to  the  interface  with  the  angle-ply  layer.  Again,  if  both  the 
90-degree  layer  and  the  free-edge  elements  were  refined,  the  improvement  of  o\  was 
not  only  on  the  90-degree  layer  but  also  on  the  entire  laminate.  Figures  (48)-(49) 
show  this  improvement. 

Comparison  of  the  interlaminar  shear  stress  t  at  the  center  line  of  90-degree 
layer  along  the  width  of  various  laminate  specimens  are  shown  in  Figures  (50M53). 
Solutions  from  both  methods  indicate  that  r  approached  finite  maximum  values  in  the 
vicinity  of  free-edge  yet  the  traction-free  boundary  condition  could  only  be  satisfied  by 
using  the  continuous  traction  0-23  element.  The  maximum  values  of  rv,  from  the 


112 


THICKNESS  Z/H 

,0.00  1 . 00  2.00  3.00  U .00  5.00  S.00  7.00  8.00  9.00  10.00  11.00 


1  Z-StRESS 
y  XZ-S1RESS 


MO.  00  60.00  80.00  100.00  120.00 

STRESSES  (PSI) 

Figure  38:  Through-Thickness  Stress  Distributions  at  the  Free-Edge  of  Speci¬ 

men  A  Rased  on  the  Continuous  Traction  Q-23  Element 


114 


-STRESS  (PSD 

25.00  50.00  75.00  100.00 


+  SPECIMEN  ft 
x  SPECIMEN  B 
0  SPECIMEN  C 
+  SPECIMEN  D 


DISTANCE  FROM  CENTER  LINE  T/B 


f-'igure  42:  Distribution  of  Z-stress  at  Mid-plane  for  Width  of  10  Plies 

Based  on  the  Continuous  Traction  Q  23  Element 


117 


S.  DO  6.00  7.00  8.00  9.00  10.00 


Z-STRESSES  (PSD 


Figure  43:  Through-Thickness  Z-stress  Distributions  of  Speimen  A  Based  on 

the  Continuous  Traction  Q-23  Clement 


118 


Through-Thic 
men  A  with 


ORE  REFINEME 
ER  REFINEMEN 


122 


ress  Distril 
•  Kefineme 


THICKNESS  Z/H 


-TZ-STRESS  (PSI) 

.00  3.00  6,00  9,00  12.00  15.00 


q  CONTINUOUS  TRACTION 
^  CONSTANT  STRAIN  ELEMENT 


126 


-TZ-STRESS 

,_p.00  5.00  10.00 

cr  i  I  _ l  ___ 


THICKNESS  Z/H 

0,00  1,00  2,00  3,00  «,00  5,00  6,00  7,00  8,00  9,00  10.00  1.1.00 


0  CONTINUOUS  TRfiCT  1  ON 
+  CONSTANT  SlflfilN  ELEMENT 


Figure  56: 


Through-Thickness  Distribution  of  YZ-stress  at  the  Center  of  the 
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Figure  57:  Through-Thickness  Distribution  of  YZ-stress  at  the  Center  of  the 

Second  Flement  From  the  Free-Hdge  for  Specimen  D 
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two  models  were  comparable.  The  constant  strain  triangular  element  approximation 
departed  significantly  from  the  Q  23  solution  in  the  vicinity  of  the  tree-edge.  This 
could  largely  be  due  to  the  nonsatisfaction  of  the  traction-free  boundary  condition. 
Figures  (54)-(57)  illustrate  through-the-thickness  stress  distribution  of  t  along  center 
line  of  the  second  element  from  the  free-edge  for  various  laminate  specimens.  The 
reason  to  choose  this  site  for  comparison  was  1  recause  of  the  singular  stress  behavior 
being  displayed  in  the  finite  element  discretization  for  both  methods  (Figures  50  5.1). 
A  close  agreement  was  generally  observed  between  these  two  solutions.  The  continuous 
traction  Q-23  analysis  show  that  t he  maximum  r  occurred  at  the  interlace  lretueen 
the  negative  angle  ply  and  the  90-degree  layers  for  all  the  specimens.  The  constant 
strain  element  does  not  have  the  capability  to  predict  this.  It  is  also  noted  that  lor 
the  same  applied  axial  stress,  specimen  D  has  tin  largest  value  of  r  ,  followed  by  A, 
B  and  C,  similar  to  the  observation  for  oq 

Table  (2)  shows  the  values  of  interlaminar  normal  stress  <rz  at  the  interface  of 
the  free-edge  for  laminate  specimens  A,  B,  C  and  D  based  on  constant  strain  element 
and  continuous  traction  Q-23  element,  respectively.  The  ratios  of  the  normal  stresses 
cry  to  the  relatively  minimal  value  among  them  are  2.11  :  1.32  :  1.0  :  2.47  for 
constant  strain  element,  and  2.12  :  1.33  :  1.0  :  2.39  for  continuous  traction  Q-23 
element.  Thus,  the  normal  stress  ratios  for  specimens  A,  B.  C  and  D,  calculated  from 
these  two  finite  element  schemes  are  comparable  with  each  other. 

Figure  (58)  shows  exaggerated  views  of  the  displacement  fields  based  on  the 
continuous  traction  Q-23  element  in  specimens  A,  B  and  C.  Figure  (59)  show's  the 
distortion  of  Specimen  D.  The  maximum  displacement  in  the  y-direction  calculated 
from  both  constant  strain  element  and  continuous  traction  Q-23  element  shown  in 


Table  2:  Comparisons  of  Delamination  Tendency  for  Various  Speimens 


(  unit  :  psi) 


Classification 

A 

B 

C 

D 

crz  from  Constant 
Strain  element 

Value 

107.74 

67.65 

51.12 

126.22 

Ratio 

2  11 

1.32 

1  00 

2.47 

a.  from  Continuous 

Stress  Element 

Value 

82  07 

bl  5m 

38.68 

92  43 

Ratio 

2  12 

1  33 

1.00 

2  39 

Note  :  Ratio= 


cr 


cr 

zmin 
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D  Laminate 


D  Laminate  Specimen 


Table  (3),  were  found  to  be  in  reasonable  agreement  for  all  these  specimens.  Also, 
from  these  figures,  it  is  observed  that  under  the  same  applied  axial  stress,  specimens  D 
had  the  largest  distortion  near  the  Tree-edge,  followed  by  laminates  A,  B  and  C  in 
descending  order.  In  other  words,  specimen  D  had  the  greatest  tendency  for  edge 
delamination  among  the  four  specimens,  and  the  delamination  would  set  in  at  the 
lowest  applied  axial  stress.  This  again  confirms  the  prediction  based  on  the 
interlaminar  normal  stress  cr  which  had  the  largest  values  for  specimen  D  as  shown 
in  Table  (2).  Based  upon  above  analysis,  we  conclude  that  both  the  distortions  and 
the  values  of  normal  stress  cr  near  the  free-edge  of  the  specimens  calculated  from  the 
continuous  traction  Q-23  element  are  consistent  with  those  of  the  constant  strain 
element  [4b]  which  is  much  more  economical  to  implement. 


1'able  3:  Comparisons  of  Maximum  Transverse  Deformation  for  Various  Specimens 


(  10  *  in) 


Classification 

Specimen 

A 

B 

C 

D 

Constant 

Strain  Element 

Value 

0.1309 

0.1032 

0.0829 

0.1935 

Ratio 

1.5789 

1.2448 

1.00 

2.3337 

Continuous 
Traction  Element 

Value 

0  1171 

0  0888 

0.0704 

0.1638 

Ratio 

1  66-,  1 

1.2619 

1  00 

2.3274 

llote  :  Ratio= - 

v 

111  in 
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6.3.1.2  Analytical-Experimental  Correlation 

An  experimental  study  was  conducted  at  AFWAL/AIFDL  [46]  to  validate  the 
analytical  results.  Part  of  the  test  results  are  summarized  in  Table  (4),  which  gives 
the  specimen  type,  strains,  and  stresses  for  the  initiation  of  edge  delamination 

determined  by  transverse  gages,  silver  ink,  acoustic  emission  and  visual  observation. 

The  table  also  contains  average  axial  stress  lor  initiation  of  edge  delamination 

calculated  front  the  continuous  traction  (J  23  element.  The  average  initial  delamination 
stresses  for  laminates  A,  11.  and  C  were  found  to  be  19.2.  23.3  and  26.4  (ksi). 
respectively  [45].  The  corresponding  finite  element  solutions  were  15.883,  25.817  and 
29.935  (ksi).  To  calculate  these  values,  the  nodal  axial  stresses  within  each  L(,'(T-12 
triangular  element  had  to  be  recovered  from  the  stresses  calculated  along  the  boundary 
of  Q-23  quadrilateral.  The  axial  loading  applied  on  each  triangular  element  could  then 
be  computed  by  integrating  the  nodal  stress  values  over  each  triangular  area.  Having 
assembled  the  element  axial  loading  over  the  whole  system,  the  average  axial  stress 
was  obtained  by  dividing  the  total  axial  loading  by  the  total  cross-sectional  area. 


Table  4:  Experimental  Results  for  Various  Laminate  Specimens 


Laminate  Specimen 

A 

B 

C 

Axial  Strain 

10  6 

6099 

3300 

2900 

Initial  Matrix  Cracking  Stress 
(  ksi  ) 

15.5 

21.7 

25.1 

Initial 
Delamination 
Stresses 
(  ksi  ) 

Gages 

20.4 

23.7 

22  1 

Silver  Ink 

17.2 

24.6 

24  9 

Acoustic 

Emission 

19.2 

23.2 

26.-4 

Visual 

18 . 1 

23.5 

25.9 

Average  Stress  Calculated 
from  Continuous  traction  Q-23 
element  (ksi) 

15.883 

25.817 

29.935 
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6.3.2  Delamination  of  1(0/-9>2/9O]s  Laminate 

A  sequence  of  tests  had  been  conducted  [51  ]  to  monitor  the  material  damage  in 
[(0/— 0)2/9O]s  laminate  specimens  under  incremental  loading.  The  value  of  0  were 

5°,  15°,  25°,  35°,  45°.  The  material  considered  here  was  T300/5208  graphite  epoxy 
with  the  following  elastic  constants: 

H|1=22xH>',psi 
H,,=  1.54xl0'psi 
G  !,=(],, =0.81  xio'psi 

*'l2  =  ,,23=()-28 

The  thickness  of  the  ten  ply  laminate  averaged  around  0.06  inch  with  width  equal  to 
1  inch.  All  these  laminate  specimens  had  also  been  analyzed  [51]  using  the  assumed 
stress  hybrid  finite  element  model  [52]  in  conjunction  with  the  quadratic  tensor 
polynomial  failure  criterion  [53]  tc  predict  the  onset  of  transverse  cracking  and 
delamination.  In  the  present  study,  continuous  traction  Q-23  element  with  a  uniform 
100  element-model  shown  in  Figure  (60)  was  used  to  analyze  one  quadrant  of  a 
typical  x=constant  plane  in  the  laminate  specimens.  The  calculated  stresses  along  the 
traction-free  edge  were  then  substituted  into  the  following  failure  criteria  to  determine 
the  possible  sites  for  initiation  of  transverse  cracks  and  delamination. 
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Figure  60: 


l(X)-Klement  Model 


6.3.2. 1  Anisotropic  Strength  and  Failure  Criteria 

With  macroscopically  homogeneous  but  orthotropic  materials,  development  of  a 
strength  theory  has  been  frequently  accomplished  by  extending  one  of  the  isotropic 
analyses  to  account  for  anisotropy.  Since  strength  theories  are  used  primarily  to 
predict  onset,  rather  than  mode  of  failure,  the  macroscopic  viewpoint  will  predominate. 
It  has  been  stated  [53]  that  all  the  failure  criteria  are  the  degenerate  cases  of  the 
tensor  polynomial  failure  criterion 

Fcr  +F  o-cr  +F ...  crcr  <x.  +....;£  1  (1941 

it  ij  i  j  ijk  i  j  k 

or,  explicitly 

r)or1+F,cr,+r3(r3+2Ft2(rlcr2+2ri3cr1cr,+21',3o-,a, 

+FllCT^F22CT^F33Cr5+F440'4+FS.S0's+F66<r6 . ^  (l95) 

Here,  or  are  the  stress  tensor  components  in  the  material  coordinates  and  F,,  F0  and  F1Jk 
etc.  are  the  components  of  strength  tensors,  all  components  are  referred  to  the  material 
principle  axes.  In  (195),  terms  associated  with  o\,  o-,,  and  tr6  which  are 
F«,  Fs.  and  F4,  are  taken  to  be  zero  since  shear  strengths  are  the  same  for  positive  and 
negative  shear  stress.  It  is  also  assumed  that  there  is  no  interaction  between  shear 
stresses  and  normal  stresses,  thus  F]6,  F?6,  F^  etc.  become  zero. 

The  strength  and  failure  criteria  considered  in  the  present  study  include  the 
maximum  stress  criterion,  maximum  strain  criterion,  Hoffman’s  criterion  and  the 
Tsai-Wu  criterion.  The  reason  for  choosing  these  criteria  was  not  only  their  popularity 


but  also  because  they  include  unequal  tensile  and  compressive  failure  strengths. 
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Maximum  stress  criterion 

Failure  of  material  is  assumed  to  occur  il  any  one  of  the  following  conditi  .s  is 
satisfied  [47] 


cr,  >X  -  cr2>VT;  ct3>Zt 
ct  ,>R;  cr.  >S;  cr  >T 

4  .1  t> 


(196) 


where  cr,,  cr,,  cr,  are  the  normal  stress  components;  cr„.  crs.  cr6  are  the  shear  stress 


components: 

V  V,. 

Z,  are  the  lamin 

a  normal  strengths 

in 

the 

respect  i\  e 

ly; 

iilld  R, 

S.  1  .are  the  si 

lear  strengths  in 

the 

vz. 

respect i\ e 

ly. 

\V  hen 

cr,.  ,cr  ..  cr,  o re 

compressixe.  they 

should 

Y(.  Y,  and 

Z, .  the 

normal  strengths 

in  compression  in 

the 

X. 

Reddy  [54]  stated  that  the  maximum  stress  criterion  could  also  be  expressed  in  the 
form  of  tensor  polynomial  criterion  as 

(cr1-XTXcr(+X<.X<r2-YTXcr2+YrXcr3-ZTXcr3+Z(.Xcr4-RX<r4+R) 


(cr— SXcr.+SXcr  — TXcr  +T)=0  (197) 

>  0  D 

(Comparing  (197)  with  (195)  and  ignoring  those  higher  order  terms,  the  strength  tensors 
are  [54] 


1  = 


1  1 


1  1 


X,  Xjr-  Y,  Y,  ’  1  3  Z, 


T  '  C 

1 


T  < 


1  J_ 
-■  ^c 


1  = - l - •  r  =  . - 1_;  F  = — 1 - 

"  XTX(  -  «  YTYr  •»  ZTZ( 

f  =  J_-  j  =_L  j-  =_L 

44  R2 '  SS  S2’  6<’  T2 
.  _  F,F2.,-  _  F,F3  .  p  _  F2F3 

f,2— y-’F.3— T’^3-  T" 
and  the  remaining  strength  constants  are  zero. 


(198) 
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Maximum  strain  criterion 

Failure  is  assumed  to  txcur  if  one  of  the  following  conditions  is  satisfied  [48] 
e,  >X  € . >  Y  •  e,>Z  T 

1  «T  2  ,1  3  .T  (]99) 

€4>r£;  65>S£:  e6>T 

where  €,,  e2,  e3  are  the  normal  tensile  strains  in  the  x,  y,  z  directions  respectively; 
e4,  es,  e6  are  the  shear  strains  in  the  yz,  xz  and  xy  planes  respectively;  Xfl,  Y(],  Z, , 
are  the  tensile  strain  strengths  in  the  x.  v,  /  directions  and  lvf.  S,.  T,  are  the  shear 
strain  strengths  in  the  vz.  xz  and  xv  planes  respectively.  Again,  expressing  the 
criterion  in  the  form  (if  tensor  polynomial  (195). 

(e  — X  X €  +  \  ,  X€.-V  ...)(€,+ Y  . .Xe .— Z  ..  Xe  +Z  Xe  — K  Xe  +R  ) 

(e.-S  Xe.+S  Xe  -T  Xe+T  )=()  (2(H)) 

>  t  >  e  t>  <  o  i 

Expressing  strains  in  terms  of  stresses  via  the  compliance  matrix  for  orthotropic 
materials,  (200)  can  be  expressed  in  the  form  of  (195)  and  we  have  [54] 

*22  °33 

2  sn  1  2  srt  1 
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Hoffman's  criterion 

Hoffman's  criterion  [49]  is  a  special  case  of  the  tensor  polynomial  criterion  for  the 
following  choice  of  the  parameters  F,  and  F,: 
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'  *  \  \  ;  1 yt  vr  ’  zT  z( 

,  =  =_1_ 

11  \r\  22  YTYC*  53  ZTZC. 
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44  R2’  55  S2’  66  T2 
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1 

1 

x,x( 

7  Z 
r< 

v,v( 

1  4- 

1 

1 

2  Z  7  V  Y  Y  \ 

^  *~i 1  t  <  vi < 


(202) 


Tsai-Wu  criteria 

The  Tsai-Wu  criterion  is  given  by 


Fcr+F  era  >1 

>  '  >j  •  j 


where 


(203) 


1  1 


1  1 


1  1 


r  i  i  p  l  j  p  i  i 

■~x7  2~y7  y^:  3"z7_z^ 
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44  R2’  M  S2’  bb  T2 
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2  V  ytycztzc 

Here,  it  is  noted  that  the  maximum  stress  and  maximum  strain  criteria  involved 
several  separate  equations,  and  there  was  no  allowance  for  interaction  of  the  stresses  or 
failure  modes.  The  Hoffman  criterion  and  Tsai-Wu  criterion,  however,  do  provide  for 
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interaction,  and  the  interaction  is  fixed.  That  is.  these  failure  expressions  are  not 


invariant  with  respect  to  coordinate  system.  As  expressed  in  terms  of  quadratic  tensor 
polynomial  shown  in  (202)  and  (204)  respectively,  the  only  difference  between  these 
two  criteria  was  on  the  determination  of  strength  tensors  F,?,  F1S,  F2J. 

63.2.2  Onset  of  Material  Damage 

luich  laminate  specimen  was  tested  individually  in  an  electro-hydraulic, 
servo  controlled  closed  loop  testing  machine  [51].  The  strain  and  nominal  stress  at  the 
first  sight  of  transverse  cracking  and  onset  of  delamination  are  summarized  in  Table 

(5). 

The  measured  strength  (ksi)  of  T300/520S  graphite  epoxv  are  given  by  [50] 

Longitudinal  tension  :  XT  =  210 

Longitudinal  compression  :  Xf.  =  200 

Transverse  tension  :  YT  =  10 

Transverse  compression  :  Y(.  =  21 

Shear  in  1-2  plane  :  S  =  13 

It  is  further  assumed  that  Zr=YT,  ZC=YC,  R=S  and  T=S/2.  Substituting  above 
information  into  (198),  (201),  (202)  and  (209)  to  calculate  F,  and  FtJ,  the  strength 
tensor  for  any  complex  stress  state  can  then  Ire  obtained  and  compared  with  the  actual 
stress  tensor.  Failure  is  assumed  to  occur  when  the  magnitude  of  the  actual  stress 
tensor  exceeds  that  of  the  strenth  tensor. 

Transverse  Cracking 

Based  on  the  stress  field  calculated  from  the  continuous  traction  Q-23  element,  the 
four  failure  criteria  discussed  in  6.3.2.1  were  applied  to  every  point  along  the 
traction-free  edge  for  all  five  laminate  specimens  at  the  strain  levels  shown  in  Table 
(5),  respectively,  when  the  first  sight  of  transverse  cracking  was  detected.  The  results 
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Table  5:  Test  Results  for  Ki0  )g  9o]^  laminate 


e 

5 

15 

25 

35 

45 

Transverse 

Strain  (*/,) 

0.493 

0.326 

0.301 

0.351 

0.532 

cracking 

Stress  (ksi) 

69.7 

39.6 

29.1 

20.5 

15.9 

Delamination 

Strain  ('/,) 

0.697 

0.383 

0.336 

0.406 

0.620 

Stress  (l:si) 

99.3 

46.7 

30.8 

23.9 

18.2 

1 49 


15.0 


63: 


Ratio  of  Stress-to  Strength  Tensor  for  [(±25)2  ■'90]^  laminate 
under  Transverse  (Tracks  Loading 


15* 


are  shown  in  Fumres  (01)(b5>  along  with  ('lion's  [50]  predictions.  Those  points  that 
lie  in  the  region  where  the  actual  stress -to-strength  tensor  ratio  is  greater  than  unity 
represent  failure.  Due  to  the  discrepancy  of  the  calculated  stress  field  based  on 
different  numerical  schemes,  the  prediction  of  the  lamina  failure  surface  from  the  same 
failure  criterion  (such  as  Tsau-Wu  theory)  varied  significantly  through  the  laminate 
free-edge.  An  obvious  failure  phenomenon  resulting  from  the  transverse  cracking 

within  the  90-degree  layer  was  detected  based  on  the  continuous  traction  Q23  element 
loi  all  the  laminate  specimens.  At  the  same  time,  initiation  of  edge  delamination  at 

the  interlaces  between  0  and  — 0  become  apparent  as  the  value  of  0  increases,  which 
was  not  indicated  according  to  Chou's  analysis.  However,  the  fact  that  transverse 
cracks  always  occurred  prior  to  delamination  in  all  cases  is  noticed,  and  this  indeed 

matches  experimental  observation  [51].  Here,  it  is  noted  that  the  magnitudes  of 

stress-to-strength  ratios  shown  in  these  figures  sometimes  departed  significantly  from 

unity  particularly  in  the  interior  of  transverse  layer  and  near  the  interfaces.  This 

could  possibly  be  due  to  the  inaccurate  insitu  transverse  strength  data  and  the 

inappropriate  assumption  of  the  interlaminar  strengths. 

Edge  Delamination 

Following  the  same  procedure  as  in  the  prediction  of  transverse  cracking,  maximum 
stress,  maximum  strain,  Hof  l  man  and  Tsai-Wu  criteria  w  ere  applied  to  every  point 

along  the  free-edge  of  various  laminate  specimens  at  the  respective  strain  level 
correponding  to  the  onset  of  delamination.  For  illustration,  failure  surfaces  predicted 
from  the  continuous  traction  Q-23  element  for  0=5°  and  0=25°  laminate  specimens  are 
showm  in  Figures  (66)  and  (67).  In  the  case  of  0=5°,  Figure  (66)  shows  that 
following  the  transverse  tracks  formed  in  the  90-degree  layer,  delaminations  were 


developed  at  the  interlaces  between  5  5  layers.  Mean  whole,  transverse  cracks  also 
extended  to  the  angle-ply  layers  under  incremental  loading.  The  fact  that  all  the 
failure  surfaces  exceeded  unity  shown  in  Figure  (67)  might  result  from  the  inaccurate 
material  strength  data.  F'or  the  0=25“  laminate,  however,  transverse  cracks  were  still 
confined  to  the  90-degree  layer  as  delamination  propagated  at  the  25/-25  interfaces. 
Based  on  these  observations,  we  can  infer  that  fiber  orientations  of  the  laminates  with 
the  same  slacking  sequence  have  placed  an  important  role  on  the  determination  of 
damage  nr sles  under  incremental  loading. 

In  general,  the  Hoffman  theorv  had  more  conservative  predictions  than  the  others 
on  the  initiation  of  transverse  cracking  within  the  90-degree  layer,  and  the  maximum 
strain  criterion  predicted  conservatively  on  the  subsequent  edge  delamination  at  the 
interfaces  between  angle-ply  laminae.  Since  the  materials  w'ere  assumed  linear  elastic 

in  the  analysis,  the  applied  strain  loading  corresponds  to  the  onset  of  transverse  cracks 
or  delamination  based  on  the  continuous  traction  Q-23  model,  and  the  failure  criteria 
would  be  expected  to  be  lower  than  the  experimental  observation.  However, 
throughout  the  analysis,  delamination  was  assumed  to  occur  in  a  state  of  generalized 
plane  strain  without  the  influence  of  transverse  cracking.  In  reality,  this  is  not  the 
case.  More  work  needs  to  be  clone  to  study  the  interrelationships  between  delamination 

and  other  damage  modes  such  as  matrix  cracking  and  filter  breakage,  etc.  Also,  many 

practical  composite  systems  actually  exhibit  extensive  nonlinear  mechanical  response  in 
shear  and  transverse  to  the  reinforcement,  resulting  in  nonlinear  laminate  mechanical 
behavior.  Extension  of  the  present  continous  traction  finite  element  procedure  to 
include  nonlinear  material  behavior,  along  with  careful  determination  of  material 

properties  and  strength  data,  may  lead  to  better  estimation  of  initiation  of  various 
damage  modes  under  incremental  loading. 
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Figure  66:  Ratio  of  Stress-to- Strength  Tensor  for  [(±5)^  /90]^  laminate 

under  Delamination  i.oading 
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SECTION  VII 


DISCUSSION 

The  problem  of  tree-edge  deiammation  in  composite  laminate  coupons  subjected  to 
unilorm  in  plane  extension  lias  been  it'.x  est  ig-.U-d.  Bel  ore  deiammation  can  be  predicted 
\  ia  a  sties'-  basts!  I  at  In  re  criterion,  an  acctiiute  siie^s  ( .dcal.oi  on  within  the  laminate, 
particularly  near  the  interlaces  .md  traction  I  ree  b.-undarv.  is  necessary.  llo\ve\er.  the 
stress  field  under  this  situation  is  highly  complex  m  nature.  Besides,  the  anisotropy 
and  heterogeneity  of  the  material  system,  and  presence  of  the  traction-free  boundary 
makes  the  analysis  difficult.  Literature  on  this  subject  was  abundant  but  an  effective 
as  well  as  reliable  solution  had  not  been  found. 

The  present  research  effort  has  resulted  in  dec  elopement  of  a  continuous  strain 
finite  element  model  in  plane  elasticity  based  on  the  compatible  cubic  interpolation 
function  proposed  by  Clough  and  Lelippa  [40].  in  which  the  normal  slope  continuity 
was  ensured  across  the  interelement  boundary.  L\ tending  the  continuous  strain  mode! 
to  the  analysis  nl  a  pseudo  t w o  dimensional  free  edge  deiammation  coupon  under 
unilorm  extension,  a  continuous  strain  held  along  both  inplane  and  transverse  directions 
was  obtained.  However,  due  to  material  anisotropy,  the  stresses  along  the  interfaces 
between  differently  oriented  layers  were  discontinuous.  Also,  traction-free  boundary 
was  not  satisfied.  The  continuous  strain  model  was  used  as  the  basis  for  the 
developement  of  a  continuous  traction  finite  element  prtxedure.  Knowing  the  fact  that 
the  displacement  field  within  each  element  iv  described  by  nodal  displacement 


components  <1111.1  their  gradients.  to  ensure  traction  lontinuirv.  a  transformation  procedure 
was  developed  to  m.ip  the  gradients  normal  to  element  Ixuindarv  to  ;i  mixed  set  ol 
degrees  ol  freedom  through  appropriate  displacement  stress  relationships.  lor  global 
assembly,  the  nodal  degrees  of  freedom  of  this  element  include  interlaminar  stress 

components  at  the  corner  nodes,  as  well  as  traction  components  at  the  mid-side  nodes 
of  each  element.  This  ensures  continuity  of  displacement  and  traction  along 
tnui  clement  boundaries  as  well  as  across  laminate  interfaces  providing  a  small 
del  ■  in. ii  ion  situation  is  considered.  At  the  same  time,  equilibrium  condition  is 
maintained  be' ween  two  ad  latent  elements  (lasers).  A  significant  aspect  of  liiis 

displacement  based  formulation  is  that  it  allows  t nation  I  ree  boundary  conditions  to  be 
specified  in  a  point  wise  sense. 

In  the  four-ply  laminate  analysis,  numerical  results  calculated  from  the  continuous 
traction  Q-23  element  generally  agreed  well  with  Pagano's  analytical  solutions  [8] 
although  these  two  schemes  were  based  upon  quite  different  theories.  For  illustration, 
'fable  (6)  outlines  the  basic  characteristics  associated  with  each  of  these  approaches. 

1  ne  approximate  solutions  for  stress  components  rx,  and  cr7,  which  play  an  important 
role  in  delannnation  of  composite  laminates,  were  calculated  using  both  approaches  and 
lound  to  have  similar  distribution.  flic  study  also  revealed  that  the  pattern  of  mesh 
rclincment  had  significant  effect  on  the  estimates  of  interlaminar  stress  field  in  the 
vicinity  ol  traction-free  edge  or  near  the  interface  between  two  differently  oriented 
layers.  Here,  it  is  essential  to  realize  that  the  continuous  traction  finite  element 
procedure  is  only  applicable  to  the  Q)-23  element  and  cannot  be  simplified  to  Q-15  and 
Q-I9  elements.  This  is  because  the  continuity  of  traction  across  laminate  interface 
cannot  be  simplified  in  the  absence  of  mid-side  nodes  at  the  interface  between  two 


Table  6:  Comparison  oT  Pay  a  no's  Theory  and  Continuous  Traction  Q-23  Element 


Method 

Pagano 

Continuous  traction  Q-23 

Variational 

principle 

Reissner’s  variational 
principle 

Minimum  potential 
energy  principle 

Type  of  formulation 

Mixed 

Displacement 

Basis  of  field 
equations 

Plate  theory 

Elastic  solid 

Assumed  inside  each 
layer  (element) 

kinematic  relations 
&  Stress  equilibrium 

Kinematic  relations 
k  constitutive  law 

Along  interlaminar 
boundary 

Continuous  traction 
&  weighted  displacement 

Continuous  displacement 
&  traction 

Assumed  stress 
v: .  r  .  t .  Z-a:-:is 

Inplane - 1 inear 

Transverse - cubic 

Quadratic 

Unknowns  in  final 
equations 

Weighted  displacements 

U  interlaminar  stresses 

Displacements 
&  interlaminar  stresses 

Solution  technique 

Direct  solving 
differential  equations 

Finite  element  method 

H.l 


adjacent  layers.  In  comparison  with  t  lie  continuous  strain  Q-23  element,  the 
introduction  of  the  transformation  process  in  the  Q  23  element  makes  the  continuous 
traction  procedure  more  expensive.  However,  the  continuous  traction  Q-23  element 
significantly  improves  the  reliability  of  the  stress  field  solution  because  of  the 
interlaminar  stress  continuity  at  the  interface  between  differently  oriented  fiber  layers, 
along  with  satisfaction  of  traction-free  boundary  condition  along  specimen  edges. 

Application  to  the  multi-ply  laminate  specimens  with  stacking  sequence  of 
I'd  —(>)..  do  ]  further  illustrated  the  potential  of  the  continuous  traction  finite  element 
procedure  in  the  ana!\sis  edge  delammalmn  problem.  Satisfactory  agreement  was 
generally  obserxed  lor  interlaminar  stress  distributions  as  well  as  laminate  displacement 
field  between  continuous  traction  Q  23  element  and  constant  strain  element  solutions 
[44]  except  that,  in  the  vicinity  of  the  free-edge,  the  constant  strain  element  was 
deficient  due  to  nonsatisfaction  of  traction-free  boundary  condition  (Ty:=0)  and  the 
assumption  of  (rxl=0)  imbedded  in  the  axisymmetric  analysis.  The  results  from  the 
continuous  traction  Q-23  element  would  be  expected  to  be  superior  to  the  constant 
strain  element  (conventional  assumed  displacement  elements)  for  prediction  of  stress 
field  in  the  free-edge  delamination  specimens.  Of  course  the  simple  axisymmetric 
model  is  economical  to  use. 

Regarding  prediction  of  damage  initiation  in  laminate  composite  coupons. 
1(0/— 6)/ 9()]s,  under  incremental  loading,  the  continuous  traction  finite  element  procedure 
along  with  some  popular  anisotropic  failure  criteria  was  found  to  be  successful  in 
modelling  some  failure  phenomena  observed  in  the  experiments.  Basically,  the  laminate 
specimens  under  analysis  were  assumed  to  be  made  of  linear  elastic  brittle  materials. 
Thus,  initiation  of  delamination  directly  led  to  catastrophic  laminate  failure  regardless 


of  the  damage  accumulation  process.  Numerical  experiment  discussed  in  the  last  section 
revealed  that  the  Hoffman  criterion  had  a  more  conservative  prediction  on  the 

initiation  of  transverse  cracking  within  the  90-degree  layer,  and  the  maximum  strain 

criterion  on  the  subsequent  edge  delamination  between  the  angle-ply  laminae  interfaces. 

In  summary,  we  conclude  that  the  proposed  continuous  traction  finite  element 
scheme  not  only  overcomes  the  drawback  of  deficient  stress  calculation  arising  in  the 
conventional  assumed  displacement  approach,  hut  also  provides  a  reliable  as  well  as 
ct  1  ec t i \ e  numerical  solution  procedure  with  a  wider  range  of  applicability  to  the 

analvsis  ol  the  1  ree  edge  delaniinaiion  problem.  Though  based  on  a  complete! \ 

different  variational  formulation,  tins  model  has  shared  the  characteristic  of  continuous 
displacement  as  well  as  traction  fields  across  laminate  interface,  with  Pagano's 

approximate  theory  derived  from  a  mixed  formulation.  Though  developed  for 

evaluation  of  stresses  in  composite  laminates,  the  continuous  traction  procedure  is  also 
applicable  to  analysis  of  layered  media  involving  material  interfaces  where  a 
two-dimensional  or  pseudo  two-dimensional  representation  is  applicable.  This  would 
include  stresses  in  layered  airfield  and  highway  pavements,  pressures  on  tunnel  lining, 


etc. 
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Appendix  A 

DERIVATION  OF  COMPATIBLE  CUBIC 
INTERPOLATION  FUNCTIONS 


This  appendix  coni. uns  a  summary  ol  lehpp.A  1 41.-12]  approach  in  «.ieri%  iny  the 
uil'ic  compatible  interpolations  for  m-plane  displacement  u  in  a  more  detailed  formal¬ 
in  order  to  derive  the  cubic  interpolation  functions  for  the  complete  triangular 
element,  three  different  coordinate  systems,  i.e.  triangular  coordinate,  local  and  global 
Cartesian  coordinates  should  be  defined  as  illustrated  in  Figure  (A.l).  The  geometry  of 
an  arbitrary  triangular  element  can  be  expressed  in  a  ('artesian  coordinate  system  by 
its  nodal  coordinates  or  its  projected  dimensions  as  shown  in  Figure  (A.2),  or 
alternately  by  its  intrinsic  dimensions  as  defined  in  Figure  (A.3). 

Lei  j  and  k  denote  the  first  and  second  cyclic  permutations  of  i=l ,2,3  (i.e.  j=2,3,l 
and  k  =  3.1.2),  the  projected  dimensions  may  Ire  written  as 

a  =  \ ,  —  x  :  h  =  v  —  V ,  (A.l) 

►  I  l  I  ■  k 

Also,  the  intrinsic  dimensions  may  be  defined  in  terms  of  the  projected  dimensions, 
deferring  to  Figure  (A.2),  define 
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/i  =  1  —  A 


( A3) 


(a  a.  +1)  b.  ) 

d  =  — — - —  ( A.4J 

I 

The  triangular  coordinates  of  any  point  ”P"  in  the  triangle  may  be  defined 

either  as  the  ratios  of  the  areas  A,  of  the  suhtnangles  subtended  by  the  point  to  the 
total  area  A  ol  the  triangle,  or  as  the  ratios  ol  the  normal  distances  n(  to  the  height 
h  .  i.e. 


A  n 

=  —  ( A  .5 ) 

A  h 

1 

as  shown  in  figure  (A.l).  It  is  noted  that  the  triangular  coordinates  are  related  by 
the  constraining  condition  £,+^+£,-1 

With  reference  to  Figure  (2),  the  displacement  interpolation  functions  for  each 
subelement  (i)  express  the  relationship  between  the  displacement  u'1’  within  the  element 
and  the  ten  displacement  components  of  its  nodal  points  r''  as  follows 

u  11  =  \<b  '  ) 1  ir '  I  1  A  .(>  I 

l  or  example,  the  nodal  displacemer.  sector  for  subelement  1  is 

Ir'V  =!u  u  ,u  ,u,,u  ,u  ,u  ,u  ,u  ,u  }  (A. 7) 

■  '2  v,  V  x,  v,  ■'  x,, 

The  corresponding  ten  cubic  interpolation  I  unctions  expressed  in  triangular  cixirdiates  are 
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where  the  subscripts  correspond  to  the  renumbered  nodes  of  the  subelement.  and  are 
the  local  triangular  coordinates  of  points  in  subtriangle  1.  With  this  convention,  the 
interpolation  functions  for  subelements  2  and  3  are  the  same  as  (A. 8)  appropriately 
permuting  the  subscripts  and  superscripts.  It  should  be  noted,  however,  that  the  nodal 
displacements  in  (A.6)  are  identified  by  node  numbers  defined  for  the  complete  element 
assembly. 

If  the  vector  {r}  of  all  nodal  displacements  of  the  complete  element  assembly  is 
written  as  the  ordered  set 


where  {<2f''|  is  similar  to  ( A. 8).  but  expanded  with  5  zeros  to  account  lor  the  nodal 
displacements  not  associated  with  subelement  1.  ami  with  appropriate  arrangement  of 
terms.  The  sectors  and  id/11}  represent  the  interpolation  functions  for  the 

external  and  internal  nodal  displacements  respectively. 

Impressing  the  displacements  in  the  other  subelements  similarly,  the  complete 
system  ol  displacements  can  be  written  as 

( A.1 1  i 

(A.ll)  is  an  expression  of  the  cubic  displacement  patterns  in  the  three  subelements. 
The  displacement  of  two  adjacent  subtriangles  are  identical  along  their  common 
boundary.  The  normal  slope  at  any  of  these  nodes  (sav  7  of  subelement  l)  is  given 
by 


1  • 
li 

I,,"' 

;  d>  '  '* 

If 

- 

:  O 

(j)  ' 

,  ,i 

u 

0.'  ' 

I  d> 

r 


(A. 12) 


where  b?',  b‘7| ’  respectively  are  values  of  { 


a#", 


at  node  7  lor  subeiement  1, 


0n  Qn 

and  n  denotes  the  axis  normal  to  the  element  boundary  To  maintain  internal  slope 
continuity,  it  is  necessary  that  01,1— — 01,1 .  where  the  negative  sign  results  from  the 
convention  that  the  positive  normal  is  directed  outwards,  lor  the  three  points  7.  8.  9. 


17-J 


The  values  of  r  which  will  satisfy  these  compatibility  conditions  are  obtained  h\ 
solving  (A.  14),  i.e. 


r  =— B  *Hr  =  Lr  (A. 15) 

t‘  O 

Substituting  the  slope  continuity  constraint  of  (A. 15)  into  (A. 11),  the  fully  compatible 
displacement  field  in  the  three  subelements  becomes 


hxplicit  expressions  of  these  functions  for  subelement  3  are 


(A.  1b) 


-  <  ,<  3-2 $ ,  sfiMK rn  ,  K ,  -M  2fj  r\  X  ,-V,C  ,1 
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<*>u  U)  = 
sl 


+  T^3^2a  l_a3^3_a2X2^  1 +  3(a  l_a3^3^2+^_3a  1  +  2a3^3  +  a2X2^3^ 


<J>u  (3)  =  { ] ( b2{ 3—b 2)+(b ,  — b3 «3  )£,£,£  3 


+  i[t’[3(-2b1+b,A(,+b:A2H|+3(b^ -b,)£2+(3b  -2b,Ml-b,X,)C3] 


«t*;3'  =  r;(3-2^1)+6A^^^1+^l3(Ml-A,)^2+(2\  -Mj)^  -3\,^] 


«>ux<3,  =  ^(ai^3-a3^)^a2-a3X3)^^2^3 


+  A^3^'^a3^3  a2^!*^ai^l"*"a3^3  ^a2^2^^a2  ai^l  ^S3^3^3^ 


<J)u.(3)  =  ^;(b^1-b1C,)+(bA-bP^l<^, 

.N  o 


+  ^ri[3(b -bfAt)^l+3(-b)Ml-b,A2+2b2)^2+(-3b2+b1p1+2b,A,X,] 


4*;  ’ ' '  -  U 3(  1  +p  X , +3(  1  +  A ,  )$  ,+<  1  -m  ,-A ,  H  J 


4< u  ' '  =  ‘ [3( a i  +  3;* . +il ,f^Ar 3( 3“ , +a 2 + , A , +( a | A . ^ ,  1 


1  7(> 


*  (-,>=i^[-3(b1+3b2+b,M^1+3(3b1+b2+b1X1)^-(b1X,-b,M2)^] 

'  3  6 

*„  <j,="<4-K«j+<K-3>] 

n4  J13 

»,  IM=-44-[r;(3{,-«,)] 

“3  I 

<t\,  ”  (A.  17) 

The  above  set  of  interpolants  is  applicable  to  all  points  lying  in  subtriangle  3.  l'or 
points  lying  in  subtriangle  1  and  2,  {<t»cl)},  {<J><2>}  can  be  written  down  by  cyclic 
permutation  of  all  subscripts  and  superscripts  in  (A. 17).  All  the  symbols  on  the  right 
side  of  (A. 17)  relate  to  the  complete  triangle. 
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LIST  OF  SYMBOLS 


A  list  of  the  most  commonly  used  symbols  and  their  general  meaning  follows. 


a,  b 


A 

A 

1 

IS 


C 


u 


s 


IJ 


s 

II 


0,, 

d 

I 

1},, 

i: 


Global  dimensions  of  a  triangle 
Area  of  triangle 
Area  of  subtriangle  i 
Material  constants 

Components  of  the  stiffness  matrix  in  the 
global  coordinate  system 

Components  of  the  stiffness  matrix  in  the 
material  coordinate  system 

Components  of  the  compliance  matrix  in  the 
global  coordinate  system 

Components  of  the  compliance  matrix  in  the 
material  coordinate  system 

Components  of  the  transformed  reduced 
stillness  maim 

Components  of  the  reduced  stillness  matrix 
Projection  of  a  corner  of  a  triangle  over  opposite  side 
Inversion  of  f) 

'i 

Components  of  the  reduced  material 
compliance  matrix 

Young's  moduli 
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i 

h 

h 

K 

11 

\ 

A 

( ' 

1) 

u" 

R 

s 

s,.  *2 

R 

n 

s 

R 

an, 

C 7 

U 

1 


('artesian  components  ol  'lie  bodv  force  vector 
Triangle  heights 

Lamina  thickness 
Llements  of  Stiffness  matrix 

Triangle  side  lengths 
Normal  to  Itoundary 

Linear  openitor  or  matrix  ol  linear  operators 
on  a  region  R 

Linear  operator  or  matrix  ol  linear  operators 
on  the  boundarv  ol  R 

Domain  ol  operator  A 

n-dinrensional  Luelidean  space 

Open  connected  region  in  1:' 

Boundary  of  R 

Complementary  subset  of  S 

Closure  of  R 

Local  cartesian  components  of  the  unit 
normal  to  a  surface 

Local  cartesian  components  of  the  unit 
tangential  to  a  surface 

Open  connected  subregion  of  R 
Boundary  of  subregion  R 
('artesian  components  of  the  stress  tensor 
Cartesian  components  ol  the  displacement  vector 
(  artesian  components  ol  the  traction  sector 


17b 


t 


s 

1 1 


0 


u 


('artesian  components  of  the  prescribed 
traction  vector 

kronecker’s  delta 

Cartesian  components  of  the  isothermal 
elasticity  tensor 

Interpolation  functions 

Assumed  displacements  in  element  m 

Assumed  generalized  nodal  displacements  at  the 
boundary  of  element  m 


Ml 

e 

k 

F 

x  y  z 

v 

G 

€ 

( 

e' 

T 

i 

e 

o 

U 

V 

w 


Assumed  generalized  nodal  displacements  internal 
to  element  m 

Vector  of  strains  in  element  m  deriving  from  u’" 
Stiffness  matrix 
Load  vectors 

Cartesian  coordinates  in  F3 
Poisson's  ratio 
Shear  moduli 

Components  of  infinitesimal  or  linear  strain  tensor 

Natural  coordinate 

Rotation  angle  from  global  axis  x 
to  material  axis  1 

Transformation  matrix 
Applied  uniform  strain  loading 

Displacement  component  in  x-direction 
Displacement  component  in  v-direction 
Displacement  component  in  /.-direction 


ISO 


u 

V 

w 

e 

fi 

n 


G 

H.  i. 
15k 


1 

'  i 


Displacement  (unction  in  x-direction 
Displacement  function  in  y-direction 
Displacement  function  in  z-direction 
fiber  orientation 
Geometric  parameters 
Linear  functional 

Displacement-stress  transformation  matrix  for  corner  node 
Displacement-stress  transformation  matrix  for  mid-side  node 
Bilinear  mapping  on  \'RXV' 

Coordinate  transformation  tensor 
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